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CHAPTER  I 


INTRODUCTION 

In  seismology  or  in  many  other  situations  signals  or  waveforms 
are  encountered  that  may  be  represented  as  the  convolution  of  a 
few  components.  The  signal  might  be  distorted  by  transmission 
through  a  linear  system.  For  example,  the  effect  of  multipath 
and  reverberation  may  be  modeled  in  terms  of  a  signal  that  is 
passed  through  a  linear  system  whose  impulse  response  is  an 
impulse  train.  One  may  wish  to  recover  the  undistorted  signal, 
or  in  the  case  of  receiving  convolved  signals  one  might  want  to 
determine  these  components  so  as  to  characterize  the  waveform  or 
physical  process  from  which  it  originated. 

The  process  of  separating  components  of  a  convolution  is 
termed  deconvolution.  In  this  process  one  must  determine  an 
appropriate  transformation  of  the  waveform  into  the  desired  com¬ 
ponent  waveforms,  [3],  [16]  . 

In  seismology  arrivals  of  various  waves  or  phases  can  be 
considered  as  the  delayed  arrival  of  more  or  less  distorted  echoes. 
The  existence  and  timing  of  these  echoes  are  not  sufficiently 
apparent  on  the  original  seismograph  traces  for  analysis.  The 
main  interest  in  this  research  is  location  and  identification  of 
target  signal  sources,  or  determining  range  from  a  target  signal 

[  source  to  a  sensor. 

: 

*  Reference  numbers  refer  to  those  entries  in  Appendix  H,  section 
2,  page  . 

L 
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Identification  is  accomplished  mainly  by  using  feature  analysis 
usually  in  the  frequency  domain;  location  consists  of  two  aspects, 
one  is  determination  of  the  angle  to  the  target,  the  other  is 
range  to  the  target.  At  the  present  time  ranging  is  accomplished 
through  use  of  doppler  techniques  or  by  using  multiple  sensors  and 
performing  cross  correlation  to  determine  the  difference  in  time 
of  arrival  of  components  of  the  signals  which  travel  at  different 
velocities  [12].  In  this  research  it  is  desired  to  use  only  one 
sensor  to  detect  both  seismic  and  acoustic  waveforms  generated  by 
the  target  signal  source  and  to  determine  the  range  by  acquiring 
the  difference  between  arrival  times  of  the  two  signals  through 
use  of  the  deconvolution  technique  -  the  cepstrum. 

One  deconvolution  technique  is  based  on  the  wiener  theory 
of  linear  filtering.  This  technique  has  been  extensively  applied 
in  processing  seismic  waveforms  in  detection  of  echoes. 

Cepstrum  analysis  ,  which  is  a  subclass  of  the  homomorphic 
deconvolution  techniques,  is  a  method  that  can  effectively  separate 
the  signals  or  determine  the  difference  in  time  arrival  of  two 
signals.  Tn  this  research  an  investigation  is  conducted  using  the 
cepstrum  analysis  technique  for  data  sources  which  are  not 
strictly  coherent  and  periodic.  The  study  has  been  divided  into 
five  parts: 

1.  The  cepstrum  is  first  defined  and  calculated  in  general. 

2.  The  cepstrum  is  then  calculated  for  noiseless  cosine 
waveforms . 

*  For  definition  of  terms  see  appendix  A. 
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3.  The  cepstrum  is  then  calculated  for  noiseless  damped  cosine 
waveforms . 

4.  The  cepstrum  is  then  calculated  for  noiseless  exponential 
waveforms. 

5.  Finally,  the  cepstrum  is  calculated  for  noisy  cosine 
waveforms . 

Cosine  waveforms  for  the  noisy  case  are  chosen  because  of 
several  small  spikes  that  usually  exist  in  the  cepstrum  plot  of 
these  type  signals,  making  them  more  prone  to  discrimination  from 
the  noise. 

The  period  of  observation  is  . 5  seconds  and  the  delay  or 
epoch  times  are  .01,  .05,  .07,  and  .12  seconds  for  the  single  echo 
case  and  .07,  .09,  and  .15  seconds  for  double  echo  case.  The 
sampling  rate  is  1/2048  seconds  and  there  are  1024  samples  in  each 
observation  interval. 

The  emphasis  is  on  a  single  echo  case  since  for  more  than  one 
echo  the  study  becomes  very  complicated  and  calculations  if  not 
impossible  are  very  cumbersome.  Furthermore  the  computer  run  time 
becomes  prohibitive  for  multiple  echoes. 

The  effect  of  the  amplitude  and  epoch  time  of  a  single  echo 
is  observed  and  in  the  noisy  case  the  study  is  done  for  several 
signal  to  noise  ratios  (SNR).  In  calculating  the  cepstrum  for 
different  cases  the  lengthy  calculations  are  left  for  the  appendix 
and  the  results  are  presented  in  the  main  text  for  the  convenience 


of  the  reader. 


CHAPTER  II 


THE  CEPSTRUM 

Definition  of  the  Problem 

In  a  listening  post  within  listening  range  of  a  target, 
it  is  possible  to  have  one  single  source  with  multipath  propagation 
or  multiple  targets  or  multiple  targets  and  multi-path.  In  either 
case,  the  received  signal  is  a  composite  of  the  individual  signals 
or  of  the  multipath  signals.  One  receives  the  signal  in  a  majority 
of  cases  in  a  form  which  is  the  convolution  of  those  generated 
signals.  The  task  may  be  to  decompose  the  signal  into  the  original 
signals.  This  process  of  decomposing  the  convolved  signal  into  the 
original  is  called  deconvolution. 

A  deconvolution  as  coimon ly  performed  by  means  of  inverse 
filtering  or  optimum  zero-lag  Weiner  filtering  suffers  from  the 
limitation  that  either  the  shape  of  the  waveform  to  be  removed 
must  be  known  or  the  assumption  that  the  wavelet  is  minimum-phase 
must  be  made  [ 32 ] . 

In  the  field  a  problem  of  prime  interest  is  one  of 
determining  only  the  range  from  a  target  to  the  sensor.  Typically, 
this  is  accomplished  by  using  two  co-located  sensors  and  performing 
an  analysis  such  as  cross  correlation  to  determine  the  difference 
in  arrival  times  between  two  signals  of  different  propagation  times 
such  as  acoustic/seismic,  EM/acoustic,  etc.  [73].  These  techniques 
require  two  sensors  and,  for  various  reasons,  it  is  desirable  to 
devise  some  technique  that  utilizes  only  one  sensor.  The  problems 
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are  numerous  in  that  all  signals  and  echoes  have  similar  properties 
and  a  technique  which  may  answer  one  problem  might  not  be  of  any 
help  in  the  other  situation. 

History 

In  1959  when  Bogert  noticed  banding  in  the  spectrograms 
of  seismic  signals,  he  realized  that  this  banding  was  caused  by 
periodic  ripples  in  the  spectra  and  this  was  characteristic  of 
the  spectra  of  any  signal  consisting  of  itself  plus  an  echo.  The 
frequency  of  these  ripples  equals  the  reciprocal  of  the  difference 
in  time  arrivals  of  the  two  waves.  Tukey  suggested  that  this 
frequency  difference  might  be  obtained  by  first  taking  the  logarithm 
of  the  spectrum,  thereby  making  the  ripples  nearly  cosinusoidal. 

In  1960,  Bogert  programmed  Tukey 's  suggestion  on  a  computer  (see 
Figure  2-1)  and  proceeded  to  analyze  numerous  earthquakes  and 
explosions.  Tukey,  noticing  similarities  between  time  series  and 
log-spectrum  series  analysis,  introduced  a  new  set  of  paraphrased 
terms.  The  spectrum  of  the  log  spectrum  was  called  "cepstrum"  [11]. 

The  first  paper  on  cepstrum  was  published  in  1963  by 
Bogert,  Tukey  and  Healy  [3].  In  this  paper  the  main  study  concerned 
seismic  waves,  but  later  on,  the  idea  of  cepstrum  analysis  was 
used  in  speech  signal  analysis  by  A.M.  Noll  [4],  in  visual  evoked 
response  (VER)  by  R.  C.  Kemerait  [31],  in  oceanography  by  J.  C. 
Hassab  [41],  and  even  in  analysis  of  full-term  and  premature 


infant's  cries  [46], 
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Analytical  Discussion 

So  far  we  are  familiar  with  the  word  "cepstrum"  but  not 
with  the  mathematical  means  to  obtain  it.  As  was  mentioned  earlier 
when  there  is  a  composite  signal  one  might  be  interested  in 
decomposing  the  signal  and  obtaining  the  original  signals  or 
simply  be  interested  in  epoch  or  delay  times.  These  two  different 
ideas  would  lead  us  in  pursuing  two  different  methods  each  suitable 
to  our  desires.  After  introduction  of  the  word  cepstrum  by 
Bogert  et_  al,  R.  W.  Schafer  [16]  introduced  another  name  obtained 
from  a  former  one  called  complex  cepstrum.  After  this  new  invention 
the  previous  word  cepstrum,  which  was  solely  used  to  show  the  method 
of  detecting  the  echo  delay  time,  was  called  power  cepstrum. 

Complex  Cepstrum 

The  complex  cepstrum  technique  was  first  described  in  a 
doctoral  dissertation  by  R.  W.  Schafer  [16]  in  November  1968.  As 
previously  mentioned  the  cepstrum  is  related  to  homomorphic  decon¬ 
volution.  The  basic  difference  is  that  a  Fourier  transform 
(magnitude  and  phase)  is  employed,  rather  than  the  power  or  energy 
spectrum.  This  is  done  because  of  the  concern  over  recovery  of  the 
actual  signals  rather  than  simply  the  different  time  of  arrivals. 

Power  Cepstrum 

As  mentioned  before  in  contrast  to  the  complex  cepstrum 
one  is  not  concerned  with  the  recovery  of  the  actual  signals  by 
decomposing  them,  one  is  interested  in  detecting  the  echoes,  and 
it  is  because  of  this  difference  that  there  are  two  different 


8 


cepstrum  techniques.  From  this  point  on  only  the  power  cepstrum 
will  be  utilized. 

Single  Echo 

In  the  simplest  echo,  values  of  time  series  y(t)  are 
multiplied  by  a  constant  a  ,  delayed  by  a  time  difference  t  and 
added  to  the  original  series  to  give  a  new  series 

z(t)  =  y(t)  +  ay (t-T )  (2-1) 

taking  the  Fourier  transform  from  both  sides  of  (2-1)  one  has 

\ 

F[z(t)]  =  F[y  (t)  +  a  y(t-x ) ] 

=  F[y(t)  ]  +  0F[y(t-i)]  (2-2) 

considering  the  time  shifting  property  of  the  Fourier  transform  one 
has 

if  g(t)  i - *  G(oj) 

then  g(t-t)  *- G(u>)  e  (2-3) 

and  then  calling  the  Fourier  transform  of  z(t)  and  y(t),  Z(w)  and 
Y(w)  respectively  one  has 

Z(u>)  =  Y(uj)(1  +  ae"jwT  )  (2-4) 

squaring  the  absolute  values  of  both  sides  in  (2-4)  would  give 

|Z(W)|2  =  |  Y  (on)  |  2  |l  +  it  e~it0T  |2  (2-5) 


taking  the  logarithm  of  both  sides  in  (2-5)  would  result 
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2  log  |Z(u)|  =  2  log  |y(oj)  |  +  log  |  l  +  ae"j“T 


(2-6) 


Parameter  a  is  the  strength  of  the  echo  signal  with  respect  to 
the  original  signal  and  it  is  usually  less  than  one.  If  the 
series  expansion  of  logg  (1  +  X)  is  examined  we  observe 


y2  y3  v4  00  (  i\n+l  „ 

loge(l+X)  =  X  -  +  iL-  _  iL.  +  ...  =  £  -  xn  _  i<x<l  (2-7) 

n=  1 


2  log  I  Z(w)|=  2  log  | Y (u)  |  +  .868590  )' 


(-1)  n  —  i  (t'T 

-  a  e  " 

n 


(2-8) 


(2-9) 


The  relation  between  logarithms  in  two  different  bases  is 


log  A 

log.  A  =  - - - 

6b  log  b 

c 


In  our  case  b  =  10  and  c.  =  e (natural  log) 


iog  A 

log10A=  logTo  =  -^295logeA 


Let  us  define  the  notation 


log  X  =  log1QX  and  In  X  =  loge  X 


then  one  can  write 


log  |  1  +  «e~j“T  |  =  .434295  [  (a  e_ja  r)n  (2-11) 


i  -1WT  I  <  i 

n  e  111 


the  complete  term  is  then 


(2-10) 


(2-12) 


V 


Taking  the  Fourier  transform  again  from  both  sides  of 


(2-12)  one  has 


Cepstrum  =  C^(z(t))  =  F[2  log  |z (w) ] ]  =  F[2  log  1 Y(w) ] ] 

+  F[. 868590  l  ane'jn“T]  (2-13) 

n-1  n 


or 

C  (z(t))  =  F [ 2  log  |Y(u>)|] 

P 

r  r  /  i\n  n  “ j na)T  -jwt  ,  (2- 

+  .868590  /  l  (-1)  a  e  e  dw 

n=l 

notice  that  here  oi  is  the  variable  of  a  function  whose  Fourier 
transform  is  desired.  Interchanging  £  with  /  gives 


C  (z(t))  = 
P 


F[2  log  |Y(u)  |  ] 


on  n+1 

+  j.  868590  l  6  ( t-n  t  )  | 

1  n 


(2-15) 


from  (2-15)  disregarding  shape  of  the  original  signal  y(t)  ,  one 
observes  a  set  of  ripples.  It  is  clear  that  the  "frequency"  of  the 
ripple  is  just  t  ,  whose  units  are  ripples  per  cycle  per  second  or 
(necessarily)  seconds  [3], 

The  log  power  spectrum  can  be  considered  as  a  frequency 
series",  if  estimated  digitally,  it  will  in  fact  be  a  discrete  fre¬ 
quency  series.  In  such  a  series  Bogert  et  al  [3],  describee 
(nearly)  cossinusoidal  ripple  bv  its  quefrcncy,  its  gamnitude,  and 
its  saphe  at  some  (frequency)  origin. 
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When  approximating  the  function  ln(l+X)  the  condition 
was  stipulated  that 

I  a  I  i  1  (2-16) 

if  a  >  1  ,  one  can  performe  some  mathematical  manipulations  as: 

i  1+  cxe  '^UT  =  l+a^  +  2  a  cos  an 

=  (1+  l/a^  +  2/a  cosut)  (2-17) 

2 

Having  a  >  1  one  can  say  the  term  1/a  +  2/a  cos  wt  for  some  t 

and  w  is  less  than  one  and  because  ln(AB)  =  In  (A)  +  In  (B) ,  one 

2 

can  filter  out  the  term  lna  and  still  have  the  ripples  t  seconds 
apart . 

As  may  be  observed  the  cepstrum  is  a  function  of  frequency 
and  is  determined  over  a  range  extending  from  zero  up  to  a  maximum 
frequency  that  is  necessarily  equal  to  the  longest  lag  in  the 
autocovariance  used  in  the  original  power-spectrum  calculations. 

The  cepstrum  is  effective  in  detecting  simple  echoes  even  when 
autocovariance  fails  completely  [3], 

In  speech  analysis  the  cepstrum  of  a  speech  signal  has 
a  peak  corresponding  to  the  fundamental  period  for  voiced  speech 
but  no  peak  for  unvoiced  speech  [A].  Cepstral  pitch  detection  has 
the  important  advantages  that  it  is  insensitive  to  phase  distortion, 
and  is  also  resistant  to  additive  noise  and  amplitude  distortion 
of  the  speech  signal.  The  cepstrum,  when  calculated  for  a  series 
of  short  time  segments  of  a  speech  signal,  has  been  shown  to  perform 
very  efficiently  as  a  pitch  determinator  particularly  for  vocoder 
applications  [ 18]  . 
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In  radar  applications  with  low  look  angles,  and  similarly 
in  sonar,  measurement  of  echo  delay  helps  reduce  interference  due 
to  reflection  by  a  neighboring  boundary.  In  the  processing  of 
bioelectronic  data,  such  as  electroencephalograms  and  visually 
evoked  responses,  the  cepstrum  helps  identify  the  origin  of 
observed  electrical  activity  [26],  [59]. 

The  non  linear  terms  ignored  in  power  series  expansion 
of  ln(l+X)  deserve  some  attention.  The  first  time  that  we  took 
the  Fourier  transform  in  (2-4)  we  had  a  term  (1+ a  e  '|a)T).  The 
absolute  value  of  this  term  when  squared  is 

l  +  2a  cosut  +  a^  (2-18) 


and  the  natural  logarithm  of  it  is 


2  2  1  2  2 
ln(l  +  2a  cos  ini  +a  )  =  (2  a  cos  wt  +  «  )  -  -  (2  a  cosiut  +  a  )  + 


2  „  2  2 
=  a  +  2  a  cos  ujt  -  2a  cos  wt 


(2-19) 


Using  trig  identities  one  has 


,  2  2  ,  2  ,1  ,  1  -  , 

2  ix  cos  uit  =  2  a  (—  +  2  cos  2  ua  ) 


’  (1  +  cos2w r) 


(2-20) 


substituting  (2-20)  in  (2-19)  gives 


2  2  2  2 

ln(l  +2  a  cos  ujt  +  rx  )  =  i  +  2  i  cos  <»u  -  a  -a  cos  2wt 


*  2  IX  cos  uit  -  >  cos  2  wi 


(2-21) 
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2 

The  second  rahmonic  represented  by  the  term  a  cos  2  ojt  contributes 
an  amount  to  the  variance  of  the  ripple  [3],  This  is  to  be  compared 

4 

with  the  contribution  of  the  fundamental.  The  ratio  is  a  /4  so 
that  rahmonics,  although  to  be  expected,  should  be  rather  small. 

Another  source  of  complexity,  would  be  provided  by  multiple 
echoes,  such  as  the  case  of  two  echoes  with  parameters  a^,  and 

a2  ’  T2  '  T^ese  would  multiply  the  original  Fourier  transform  by 

1  +  a  i  e"j“Tl  +  a2  e_jwT2  (2-22) 

and  the  power  spectrum  by 

2  2 

1 +2a^  cosojt^  +  2  a2  cos  a)T2  +  2a^a2  cos  io(Tj- T2) +a^  +  a2  (2-23) 

so  that  the  amount  added  to  the  log  spectrum  is 

2  2 

2  cos  on  ^  +  2a2  cos  un2  -  a  ^  cos  2  an  «2  cos  2  un2 

-  2a^a  2cos  oj(t^  +  t (2-24) 

to  quadrature  accuracy.  Thus  in  addition  to  ripples  at  and  t 2  , 
one  expects  a  ripple  at  frequency  (t^  +  >2)  with  gaminitude  of 
order  twice  that  of  the  rahmonic  at  quefrencies  2\  and  2t2  . 

In  (2-24)  if  higher  order  approximation  is  used,  the  difference 
ripple  would  show  up.  Appendix  B  shows  complete  work  for  multiple 


echo  case. 


CHAPTER  III 


THE  SHORT  TIME  AVERAGE  CEPSTRUM 

The  terminology  'short  time  average  cepstrum'  might  be 
confusing  or  in  some  cases  misleading,  but  if  the  reader  bears  with 
us  for  a  short  time  the  ambiguities  will  be  clarified.  In  the 
previous  chapter  the  cepstrum  was  discussed  in  general,  and  the 
period  of  integration  for  the  Fourier  transform  was  from  -°°  to  . 
Most  of  the  realizable  experiments  are  causal,  that  is,  they  cannot 
exist  for  the  time  period  to  0.  The  use  of  the  digital  computer 
for  analysis  adds  another  restriction,  limiting  the  length  of  the 
experiment  even  if  it  is  performed  in  real  time. 

In  most  applications  one  is  interested  in  observing  an 
event  in  a  limited  time,  or  in  the  other  words,  observing  the  event 
through  a  window.  If  the  cepstrum  is  to  be  evaluated  for  a  composite 
signal  in  the  aforementioned  circumstances;  the  limits  of  integration 
must  be  changed,  and  it  is  because  of  these  limitations  that  the 
result  is  called  the  short  time  averaged  cepstrum. 

To  derive  some  analytical  expressions  for  short  time 
average  cepstrums  we  choose  some  functions  of  interest  and  find 
their  cepstrum.  As  mentioned  in  the  introduction  three  types  of 
signals  have  been  chosen: 

1.  cosine  signal 

2.  damped  cosine  signal 

3.  damped  exponential  signal. 

The  cepstrum  expression  is  derived  for  the  cosine  and  the  damped 


exponential  signals  only  when  we  have  single  echoes.  The  case  of 
damped  cosine  is  a  special  form  of  1.  or  3.  depending  on  the 
damping  coefficient. 
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Cosine  Signal 

The  choice  of  signal  has  been 


g(t)  =  cos  a)  (t-T-)u(t-T0)+a  cos  0)  (t  -  T,)u(t-T.) 

o  1  2.  o  3  4 

(3-1) 


which  is  shown  in  Figure  3-1.  As  one  observes  from  the  figure  the 
signals  are  not  in  phase  and  they  begin  at  two  different  times. 

The  complete  calculation  is  left  for  Appendix  C,  but  some  of  the 
results  are  mentioned  here. 

Taking  the  Fourier  transform  from  both  sides  of  (3-1) 


gives 

F[g(t)] 


F[cos  0)  (t-T1)u(t-T.)  +  a  COS  U)  (t-T  )u(t-T.)] 

O  1  l  O  j  4 

T  .  T  . 

/  coswQ(t  -t -  )  e  ^Wtdt  +  j  cos  wo  (t  -t  )e  dt 
t  2  "t4 

(3-2) 


taking  logarithm  of  the  spectrum  found  from  this  Fourier  transforma¬ 
tion  and  after  some  mathematical  manipulations  one  encounters  terms 
with  the  following  forms 


cos  i,!  (t  -  r^)  and  cosw(t-t^). 

Terms  containing  cosine,  would  have  impulses  in  their  Fourier 
transforms,  so  in  our  case  presence  of  ripples  at  t  -  and  r  -  i 
is  obvious.  Because  of  the  symmetrical  property  of  the  Fourier 
transform  one  should  also  obtain  ripples  at  ^  and  ,  and 
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at  T  "  t2  an<*  T  ~  t4  •  Figure  3-2  we  show  such  ripples  and 

also  existing  ripples  of  the  harmonics  of  t  ^  and  x^  . 

Throughout  the  calculation  the  window  or  time  frame  was 

from  t  =  0.  through  t  =  x  seconds,  and  the  peaks  are  at  t ^  and 

x.  .  Now  if  the  time  frame  is  T2  -  T1  seconds  where  t,  <  x„<  T1 
4  12 

and  x^  <  T2  <  x  as  shown  in  Figure  3-3,  we  show  by  analysis  in 
Appendix  D  and  observe  in  Figure  3-4  that  the  only  peak  is  at 
x  =  T2  -  T1  or  its  multiples. 


Damped  Exponential  Signal 

For  the  damped  exponential  case  the  following  function  is 

used . 

g(t)  =  (t-x1)e  C^t -Tl^u(t  -  x^  +  a(t  -x  2>e  C^t  T2^u  (t  -x2>  (3-3) 

which  is  shown  in  Figure  3-5.  Taking  Fourier  transform  of  both 
sides  of  (3-3)  gives 

x 

F  [  g  ( t )  ]  =  /\t- x  )  e"c(t"Tl)  e'jwt+  (  (t -t  )  e~c(t"T2)e_ja,tdt 

T1  2 

(3-4) 

when,  after  some  mathematical  manipulations,  the  Fourier  transform 


*Keys  for  the  figures  are 

1  WNDW  =  1  windowing  is  preformed 

I  WNDW  =0  no  windowing  is  performed 

’  WNDF  =  1  smoothing  is  performed 

T  WNDF  =0  no  smoothing  is  performed 

-  - ^  _  'I’otal  number  of  sample  points  in  the 

original  composite  signal 


m -  -  '  — aJUiiHlil 


.*cB3£i-D3«H«tr-raui!  i  Kim-rmni  vs.  #cas tsoo^pi * tr-moaj  i  stm-musi 
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Figure  3- A.  QUEFRENCT  SECONDS  t.  =  .07  T?  =  .12 


K 


S 


0 


TIME  SECONDS 

IT-THUll  EXF  r-SOa-TfUUI  I  UCr-THUll  *.  6  (f- IRL*2I  EXF  (-30  t T- rHPct  I  im-THUc’l 

Figure  3-5.  a)  A  damped  exponential  waveform, 

h)  its  echo,  c)  the  composite  signal 


i  *u  cr-ram  i  ♦.  Ex  't-rauai  *e:<p  c-sq.  * (r-  :al'2(  i  *u  cr-rauat 


QUEFRENC1T  SECONDS 
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has  been  calculated  and  the  magnitude  of  the  power  spectrum  is 
found  one  encounters  terms  such  as 

A  cos  0)  (  T  -  t ^  )  ,  B  sin  w  (  t  —  t^)  ,  C  cos  u  (  t  -  t^)  , 

D  sin  a)  (  t  -  T£)  >  E  cos  oi(  ~  T^)  »  and  F  sin  w(x2  -  x^) 

Taking  again  the  Fourier  transform  of  the  log  spectrum  of  (3-4), 
one  has  to  take  Fourier  transform  of  a  few  terms  containing  cosine 
and  sine  waveforms  of  and  T]_  »  which  would  result  peaks 

at  x^  ,  x 2  and  x^  -  x^  .  But  as  is  shown  in  Appendix  E,  the 
magnitude  of  the  ripples  at  x  and  T2  is  very  small  and  in  the 
cepstrum  plot  of  Figure  3-6  it  is  shown  that  it  is  not  recognizable 
(it  is  safe  to  say  that  with  the  choice  of  c  and  a  one  might  have 
different  degree  of  distinguishability  at  x^  and  t2>  .  In  fact  we 
show  in  Appendix  E  for  choice  of  c = 30. ,  t  =  .5,  =  .07  and 

T2  =  .12  seconds  and  a  =  .6  this  relative  ratio  is  of  the  order 


Damped  Cosine  Signal 

In  the  damped  cosine  case  the  equation  to  be  used  is 

g(t)  =  (t  -t  )  cos  oj  (t  -  t  )e  C^t  Tl^  +  ( t  —  t  )  cos  u)  (t-T  )eC^t  T  2^ 

1  o  1  2  o  2 

(3-3) 

which  is  shown  in  Figure  3-7.  Analytical  calculation  of  the  cep¬ 
strum  for  this  case  is  extremely  difficult,  but  experimental  results 
have  shown  the  cepstrum  behaves  like  the  cosine  or  damped  exponential 
case,  depending  on  values  of  t,  and  x ^  .  In  the  next  chapter 


the  experimental  results  are  discussed. 


In  all  three  cases  If  one  takes  into  account  the  higher 
order  terms  in  ln(l  +  X)  expansion,  ripples  are  observed  at 
+  t 2  and  its  harmonics.  Appendix  C  shows  the  existence  of 
such  ripples,  at  the  points  and  T2  +  Ti  *  Figures  3-7 

and  3-8  show  damped  cosine  waveform  and  its  cepstrum  respectively. 


Cepstrum  for  a  DAMPED  COSINE  Waveform  With 


CHAPTER  IV 


EXPERIMENTAL  RESULTS 

As  previously  stated  the  objective  is  to  investigate  the 
feasibility  of  applying  cepstral  analysis  to  echo  or  epoch  detec¬ 
tion.  The  study  is  conducted  for  two  cases  of  noiseless  and 
noisy  conditions,  with  single  or  multiple  echoes.  The  computer 
study  is  performed  using  Univac  1108  computer  at  Mississippi  State 
University  Computing  Center. 

Noiseless  Case  (Single  Echo) 

For  the  study  of  single  echo,  three  types  of  signals; 
cosine,  damped  exponential,  and  damped  cosine  waveforms  were  chosen. 
The  choice  for  cosine  waveform  is 

cos(2007r(t  -  x^DuCt  -  t^)  +  .6  cos  (200ir(t  -  T2))u(t  -  t  )  ^  ^ 

The  above  equation  is  comprised  of  two  parts,  the  first  part  being 
the  original  signal  and  the  second  part  the  echo.  If  is  not 

zero  one  has  the  case  of  epoch  detection  for  the  original  signal 
and  echo  or  eopch  detection  for  the  second  signal. 

The  total  length  of  the  time  window  is  .5  seconds.  The 
signals  are  sampled  at  2048  samples  per  second  and  the  delay  times 
are  =  0.0,  =  .01,  and  .05  seconds.  Also  the  cases  of 

t  =  .07,  i „  =  .12,  and  t.  =  .21  seconds  are  discussed.  For  each 

1  i.  I 

case  of  t ^  and  T2  the  study  is  conducted  in  four  different  manners: 

1.  with  no  windowing  or  smoothing  (NW,NS) 

2.  with  no  windowing  but  smoothing  (NW,S) 
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3.  with  windowing  but  not  smoothing  (W,NS) 

4.  with  windowing  and  smoothing  (W,S) 

In  case  1,  when  no  windowing  or  smoothing  is  used  the 
signal  is  sampled,  its  Fourier  transform  and  consequently  its 
power  spectrum  is  found  and  then  the  Fourier  transform  of  logarithm 
of  this  power  spectrum  is  calculated  to  give  the  new  power  spectrum 
or  the  cepstrum.  In  case  2,  the  windowing  is  performed  in  the 
time  domain,  before  the  first  Fourier  transform  is  calculated.  The 
windowing  process  is  based  on  Kaiser-Bessel  equations  with  8=8. 

In  case  3,  the  smoothing  process  only  is  used.  The  smoothing 
process  is  based  on  the  Hanning  method.  The  smoothing  is  performed 
after  finding  the  log  spectrum  and  when  ready  to  start  the  Fourier 
transform  operation  again.  In  case  4,  both  windowing  and  smoothing 
are  done. 

The  results  are  illustrated  in  Appendix  G  showing  that 

neither  windowing  nor  smoothing  was  of  any  help  in  detecting  the 

echo,  or  epoch  delay  time.  Tables  4-la  and  4-lb  show  what 

effects  the  use  of  windowing  and/or  smoothing  create.  In  other 

parts  of  the  study  x^  was  chosen  to  be  a  multiple  of  so 

again  is  one  of  the  harmonics  of  .  The  result  is 

obvious,  as  one  sees  a  higher  peak  is  at  only,  with  no 

significant  difference  in  the  value  of  the  cepstrum  at  other  points. 

These  results  are  shown  at  Table  4-2  and  can  be  compared  with  the 

values  at  Table  4-lb.  The  reason  for  this  increase  in  value  of  the 

cepstrum  is  that  instead  of  having  only  the  amount  of  cepstrum  at 

this  point,  harmonics  of  the  cepstrum  at  point  add  to  the 

original  cepstrum. 
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Included  in  the  study  was  variation  of  the  cepstrum  as 
a  or  x 2  or  both  were  changed.  Tables  4-3a  and  4- 3b  show  variation 
of  cepstrum  for  a  fixed  x ^  and  variable  a  .  As  may  be  seen  the 
value  of  the  cepstrum  remains  almost  the  same  in  the  range  of  a  =.2 
to  a  =  1.6.  Table  4-3b  shows  the  normalized  values  of  cepstrum  in 
Table  4-3a.  In  table  4-4a  the  value  of  the  cepstrum  at  some  impor¬ 
tant  points  is  presented  and  as  may  be  seen,  changing  a  changes 
the  value  of  the  cepstrum  (notice  here  is  no  more  zero). 

These  results  comply  with  what  was  found  in  Appendix  C.  Table  4-4b 

k 

is  the  normalized  form  of  Table  4-4a  . 


TABLE  4- la 


POWER 

CEPSTRUM 

AT /FOR 

NW,NS 

NW,S 

W,NS 

W,S 

- - 

N5 

II 

o 

h-» 

381.339 

339.922 

54.4870 

6.3023 

t2  =  •°5 

415.792 

290.384 

47.9950 

4.9561 

Effect  of  using  windowing  and/or  smoothing  on  the  power  cepstrum 
of  cosine  waveform,  see  Figures  in  Appendix  G  -la,  2a  for  time  wave¬ 
form  and  Figures  in  Appendix  C-lb,  2b  for  cepstrum  plots.  shows 

the  location  of  the  highest  peak  value  as  well  as  being  the  epoch  or 
delay  time  for  the  echo  signal.  =  0  for  both  cases. 


*  Key  to  the  tables 

NV:  means  no  windowing  is  performed 

W:  means  windowing  is  performed 

NS:  means  no  smoothing  is  performed 
S:  means  smoothing  is  performed 


POWER 

CEPSTRUM 

AT/FOR 


NW,NS 


NW,S 


W,NS 


T1 

265.834 

38.6806 

6.3689 

.5999 

T2 

154.062 

31.1605 

5.2801 

.3493 

BB 

24.409 

2.61198 

T2  +  T1 

37.1550 

5.5992 

g 

Value  of  the  power  cepstrum  at  the  epoch  times  as  well  as  at 
the  difference  and  sume  of  the  epoch  times  for  =  .07  and 
=  0.12  seconds  for  a  cosine  type  waveform  (see  Figure  G-3), 
also  of  notice  is  the  effect  of  windowing  and/or  smoothing. 

TABLE  4-2 


POWER 

CEPSTRUM 

AT/FOR 

NW,NS 

NW,S 

W,  NS 

W,S 

ri 

264-2950 

9.6876 

1  .  3447 

.1123 

'2 

158-4480 

18.6844 

3.4655 

.  1072 

T  2  "  T  1 

76.9708 

27. 1098 

‘ 

4. 1172 

.  2226 

rl+  T  2 

95.4120 

41.1492 

7.4753 

.1962 

This  table  is  similar  to  Table  4-lb  except  i  =  .21  seconds 
a  multiple  integer  of  t  =  .07  seconds.  See  Figures  C-4a  and 
C-4b  for  the  waveform  and  its  cepstral  plot. 


TABLE  4- 3a 
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The  damped  exponential  signal  with  a  single  echo  is 

(t-  T1)e'30(t‘Tl)u(t-T1)+  .6  (  t  -  t2)  e_30(t_T2)u(t  -  x2)  (4-2) 

with  equation  (4-2)  being  comprised  of  two  parts,  the  first  part 
the  original  signal  and  the  second  part  the  echo.  The  damping 
coefficient  has  been  chosen  so  as  to  drive  the  signal  down  in  a 
way  which  approximates  some  physical  phenomenon  which  we  are 
modeling.  As  in  the  previous  case,  the  study  is  done  for  four 
different  conditions  depending  on  the  use  of  windowing  and/or 
smoothing . 

Table  4-5a  shows  variation  of  the  power  cepstrum  and  effect 
of  windowing  and/or  smoothing  on  the  power  cepstra  when  =  0 
and  t2  assumes  two  different  values.  Table  4-5b  shows  the  power 
cepstrum  at  four  important  points  when  i  =  .07  and  t  n  =  .12 
seconds  with  a  =  .6.  Table  4-5c  shows  the  peak  value  at  the 
same  four  locations  as  for  Table  4-5b  but  for  this  table  x?  is 
an  integer  multiple  of  the  t  .  As  mentioned  in  the  case  of  cosine 
signals  windowing  or  smoothing  would  not  do  any  good,  but  rather 
it  actually  degrades  the  results. 

Looking  closely  in  Table  4-5b  one  rinds  it  interesting 

to  see  that  the  highest  peak  value  of  the  power  cepstrum  is  located 

at  T2  _  rf  ’  t^’e  same  Placc  as  was  mentioned  in  Chapter  Ill.  for 

this  Table  C,  the  damping  coefficient,  is  equal  to  30.0.  As  was 

previously  noted  the  location  of  the  highest  peak  value  is  very  much 

related  to  the  damping  coefficient.  In  Appendix  C  it  may  be  seen  how 

changing  C,  the  damping  coefficient,  changes  the  local  ion  of  the 
highest  peak  or  its  value  in  the  cepstrum  plot. 


TABLE  4-5a 


POWER 

CEPSTRUM  NW,NS 
AT /FOR 


NW,S  W,NS 


t ^  =  .01  173.026  41.6319  I  5.3620 


x 2  =  .05  63.347  20.6159  !  2.3575 


Effect  of  using  windowing  and/or  smoothing  onthe  power  cepstrum 
of  damped  exponential  waveform,  shows  the  location  of  the 

heighest  peak  value  as  well  as  being  the  epoch  or  delay  time  for 
echo  signal.  =  0  for  both  cases.  For  the  waveforms  and  their 

cepstral  plots  see  Figure  G-ll  and  G-12. 


TABLE  4- 5b 


POWER 

CEPSTRUM 

AT/FOR 

NW,NS 

NW,S 

T1 

.5413 

3.7263 

1 2 

2.3317 

. 

.  7425 

T  2~  1 1 

257.817 

28.2807 

T2  +  rl 

.  2384 

.4015 

_ 

128]  .00981 


Value  of  the  power  cepstrum  at  the  epoch  times  as  well  as  at 
the  difference  and  sum  of  the  epoch  times  for  :  ^  =  .07  and 
x ^  =  0.12  seconds  for  a  damned  exponential  waveform,  also  of 
notice  is  the  effect  of  windowing  and/or  smoothing.  See 
Figi  res G-l 3a ,  13b,c,d,e,  for  the  waveform  and  its  cepstral 
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=  0  and  t 2  assumes  two  different  values.  Table  4-6b  shows 
the  power  cepstrum  at  four  important  points  when  =  .07  and 
r ^  =  .12  seconds  with  a  =  .6  .  Table  4-6c  shows  the  peak  value 
at  the  same  four  locations  as  for  Table  4-6b  but  for  this  table 
t 2  is  an  integer  multiple  of  .  As  was  mentioned  for  the  two 

previous  cases,  windowing  or  smoothing  would  not  do  any  good,  but 
rather  it  actually  degrades  the  results. 

Looking  closely  at  table  4-6b  we  notice  that  the  behavior 
of  this  case  is  like  the  behavior  of  two  previous  cases  combined 
together  meaning  that  we  have  noticeable  peaks  at  t  ,  r2  and  r2  -  t  . 
For  this  Table  c,  the  damping  coefficient  is  equal  to  2  and  the 
frequency  of  cosine  signal  is  100  Hz  same  as  for  cosine  type  case. 

As  we  said  earlier  location  and  the  height  (strength)  of  the 
highest  peak  very  much  is  dependent  on  the  damping  coefficient. 

In  Appendix  G  we  will  see  how  changing  c,  the  damping  coefficient, 
has  changed  the  shape  of  the  cepstrum  plot. 


TABLE  4 -6a 


POWER 

CEPSTRUM 

AT /FOR 

NW,NS 

NW,S 

r  | 

I 

W.NS  | 

1 

i 

W.S  i 

r  2  =  .01 

326.876 

177.362 

: 

28.9025  1 

3 '<9  3  - 

r2  *  .05 

292.795 

196.807 

3  3,611.'’  ) 

'em  3 

Effect  of  using  windowing  and/or  smoothing  on  the  power  eepstrnni 
of  damped  cosine  waveform.  r.  shows  the  location  ot  the  In-iglnst 
peak  value  as  well  as  being  the  epoch  or  delnv  tine  lor  n  lie  signal. 
Tj  =  0  for  all  two  cases.  For  the  waveforms  ami  their  cepstral 
plots,  see  Figures  0-6  and  G-7. 
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Multiple  Echoes 

Because  of  complexity  In  the  case  of  multiple  echoes  we 
restricted  our  study  for  two  and  three  echo  cases.  As  we  can  see 
from  Figures  4-1  through  4-4,  when  the  number  of  echoes  increases, 
it  becomes  very  difficult  to  detect  the  delay  time  of  these  echoes. 
In  study  of  multiple  echoes  for  the  cosine  waveforms  we  tried  to 
change  the  time  between  two  echoes  and  see  the  effect  on  the  power 
cepstrum.  As  we  see  in  Table  4-7a,  if  a  Is  fixed  there  is  no 
change  in  the  magnitude  of  the  power  cepstrum,  at  the  two  epoch 
times,  but  the  value  of  the  cepstrum  changes  at  the  sum  and  the 
difference  of  the  epoch  times.  This  will  further  show  the 
dependence  of  the  power  cepstrum  on  epoch  times  as  well  as  other 
parameters.  Table  4-7b  shows  the  normalized  values  of  Table  4-7a. 
For  more  plots  of  multiple  echo  case  see  Appendix  G. 

Table  4-8,  4-9  and  4-10  show  the  magnitude  of  the  cepstrum 
at  some  important  points  for  two  echo  cases  also  reflects  the  effect 
of  windowing  and  smooting,  notice  we  have  peaks  at  the  time 
difference  as  well  as  sum  of  delay  times.  This  is  exactly  what 
we  expect  from  a  cepstrum  plot. 


Cepstrum  Plot  for  Single  Echo  Cosine  Waveform 


Cepstrum  Plot  for  Single  Echo  Damped  Cosine  Waveform 


Plot  for  Single  Echo  Cosine  Waveform 


0,13  0,2?  0,31  0,38 

QUEFRENCY  SECONDS  T  =  .07  T  =  .12 


Fioure  4- 3a .  GUEFRENCT  SECONDS  x  =  .07  t  =  .09  t  =  .15 
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Three  Echoes  Plus  the  Original  Signal 


TABLE 
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TABLE  4-8a 


CEPSTRUM 

AT /FOR 

NW,NS 

NW,S 

— 

W,NS 

W,S 

T1 

228.4860 

103.2690 

17.2839 

1.6130 

T  2 

131.9120 

44.8394 

7.0003 

.4590 

mm 

68.2902 

15.5707 

1.4308 

.0430 

14.5433 

113.188 

20.6583 

2.1491 

HB 

131.9120 

44.8394 

7.0003 

.4590 

t3  '  t2 

228.4860 

103.2690 

17.2839 

1.6130 

BBS 

68.2902 

1.4308 

.0430 

T1  +  T  3 

25.2834 

1.0093 

.  5299 

.0480 

'2  +  T  3 

68.2902 

15.5707 

1.4308 

.04  30 

Tl  +  T2+  T  3 

131.9120 

44.8394 

7.0003 

.4590 

Cepstrum  values  at  epoch  times,  sum  and  the  difference  of  epoch 
times  for  a  signal  comprised  of  three  cosine  waveforms.  The  epoch 
times  are  =  .07,  1 0  =  .12  and  ' 3  =  *19  seconds.  Notice  the 
difference  between  two  epoch  times  is  another  one.  Also  notice 
the  effect  of  windowing  and/or  smoothing. 


* 

1 


CEPSTRUM 
AT /FOR 


NW,NS 


T1 

25.9059 

4.8733 

.6452 

.0566 

T  2 

110.4580 

7.4845 

.8232 

.0527 

T  3 

13.4191 

.4786 

.0504 

.0058 

T2  "  T1 

223.8570 

25.5724 

2.9127 

.2899 

T3  "  T1 

110.4580 

7.4845 

.8232 

.0527 

7  3  '  7  2 

223.8570 

25.5724 

2.9127 

.2899 

T1  +  T2 

13.4191 

.4786 

.0058 

T1  +  T  3 

12.7848 

.5522 

.0707 

.0057 

1 2  +  r3 

13.4191 

.4786 

.0504 

.0058 

11+  l2+  r  3 

2.2451 

.3971 

.0614 

.0080 

Cepstrum  values  at  epoch  times,  sun;  and  difference  of  epoch  times 
for  a  signal  comprised  of  three  damped  exponential  waveforms.  The 
epoch  times  are  x^  =  .07,  =  .12  and  ;  =  .19  seconds.  Notice 

the  difference  between  two  epoch  t imes  is  another  one.  Also  notice 
the  effect  of  windowing  and/or  smoothing. 


TABLE  4- 10a 


Cepstrum  values  at  epoch  times,  sum  and  the  difference  of  epoch 
times  for  a  signal  comprised  of  three  damped  cosine  waveforms.  The 
epoch  times  arc  =  .07,  =  .12  and  t  ^  =  .19  seconds.  Notice 

the  difference  between  two  epoch  times  is  another  one.  Also  notice 
the  effect  of  windowing  and/or  smoothing. 


t2 

77.9314 

T3 

36.1534 

T2-  T 

■ 

178.248 

t3  ~  T: 

■ 

65.9162 

t3  "  T 

2 

30.5698 

42.2664 
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Noisy  Case 

The  study  included  the  cosine  type  signals  for  the 
noisy  case.  The  reason  for  using  the  cosine  waves  was  that  in  the 
cepstrum  plot  of  a  cosine  waveform  there  are  many  ripples  besides 
the  original  ones,  and  this  makes  a  good  choice  for  a  noisy  case. 
Different  signals  to  noise  ratios  were  chosen  from  2.43  up  to  15.43 
db  and  it  was  found  that  the  cepstrum  plot  is  not  of  any  help  up  to 
15  db.  It  is  for  signal  to  noise  ratios  above  this  level  that  the 
epoch  or  delay  times  in  the  cepstrum  plots  are  distinguishable. 

Figures  4-5  and  4-6  show  the  cepstrum  plot  for  signal  to 
noise  ratios  of  2.43  and  15.44  db.  As  in  all  of  the  previous  cases 
there  is  no  discernable  improvement  in  the  cepstrum  plots  by 
using  windowing  and/or  smoothing. 
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CHAPTER  V 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  power  ccpstrum  technique  demonstrated  the  ability  to 
determine  the  epoch  times  of  the  echoes  for  the  single  and  multiple 
echo  cases  with  echo  amplitude  being  less  or  greater  than  that  of 
the  reference  wavelet.  It  was  noticed  that  changing  the  relative 
magnitude  of  the  echo  with  respect  to  the  original  signal  would  not 
change  significantly  the  peak  magnitude  of  the  power  cepstrum.  The 
experimental  results  proved  to  be  in  compliance  with  the  theory. 

It  proved  once  more  as  reported  in  the  literature,  that  the  cepstrum 
is  very  much  data  related  as  seen  in  the  cases  of  damped  exponential 
and  damped  cosine  waveforms.  In  the  damped  exponential  type  signals 
it  war  apparent  that  if  the  relative  magnitude  of  the  signal  is 
small  compared  to  the  original  one  and  also  if  the  damping  coeffi¬ 
cient  is  high  enough,  the  only  place  that  a  strong  peak  occurs  in 
the  cepstrum  plot  is  at  the  time  difference  of  the  echoes,  while 
in  the  case  of  the  cosine  type  waveforms,  only  at  the  epocli  times 
is  the  cepstrum  peak  strong  enough  to  be  detected  and  the  time 
difference  is  not  discernable.  In  the  case  of  damped  cosine  a 
shift  was  noticed  in  the  form  of  the  power  cepstrum.  If  the 
damping  coefficient  is  high  it  will  react  like*  a  damped  exponential 
waveform,  and  when  damping  coefficient  is  not  high,  it  will 
behave  as  a  cosine  waveform. 

When  both  signal  and  its  echo  start  at  anv  point  in  time  except 
at  the  time  zero,  and  i  is  changing,  the  value  of  the  cepstrum  at 
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two  different  epoch  times  would  change,  but  sum  of  these  peaks 
together  will  remain  almost  constant  (see  Table  4-4b  for  cosine 
type  signal). 

Adding  some  noise  to  the  signals,  the  cepstrum  was  found  to 
be  detectable,  although  the  signal  to  noise  ratio  (SNR)  should  be 
greater  than  15.0  db  (only  the  case  of  cosine  waveform  was  tested). 

Complexity  and  enormous  number  of  peaks  at  the  epoch  times 
and  their  harmonics  makes  it  practically  impossible  to  deal  with 
more  than  two  echoes,  as  seen  in  our  plots  for  three  echoes,  even 
in  the  noiseless  case  it  is  very  hard  to  exactly  detect  the  correct 
epoch  times  of  the  echoes.  In  double  echo  case,  no  detectable 
change  in  the  magnitude  of  the  peaks  in  the  power  cepstrum  were 
noticed  when  the  time  distance  was  kept  between  two  echoes  constant. 

The  study  was  conducted  with  computer  simulation  techniques, 
which  might  be  different  from  real  time  problems,  but  comparison 
with  the  theoretical  results  showed  that  it  is  not  very  far  from 
being  correct.  Most  of  the  previous  studies  which  have  been  made 
using  cepstrum  techniques  have  been  in  the  field  of  seismic  or 
oceanography  and  little  has  been  done  in  medicine  with  the  excep¬ 
tion  of  a  few  theses  and  dissertation  researcii  from  the  University 
of  Florida.  Further  study  should  be  conducted  in  the  neurological 
field,  where  reflection  plays  an  important  role.  If  the  technique 
is  improved  enough  and  progress  is  made,  then  one  day  human  beings 
might  suffer  less  by  early  detection  of  any  malfunction  in  their 


neurological  system. 


In  the  battlefield  if  the  study  is  conducted  for  a  majority 
of  armored  vehicles,  perhaps  by  making  a  table  look-up  type  device, 


the  type  of  the  vehicle  can  be  recognized  and  so  its  capture  or 
destruction  would  be  possible.  Today  with  the  help  of  LSI  technol¬ 
ogy,  it  is  possible  to  make  FFT  processors  and  so  it  is  relatively 
easy  to  build  processing  equipments  in  which  cepstrum  techniques 
can  be  of  help. 

Finally  in  Tables  5-1,  5-2,  and  5-3,  is  a  quick  cross 
reference  for  the  study. 


* 


k. 


TABLE  5-2  (MULTIPLE  ECHOES) 
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DEFINITION  OF  TERMS  USED  IN  THIS  PAPER 


In  order  to  get  a  meaningful  terminology  for  the  cepstrum 
technique,  Bogert  e_t  aQ  suggested  the  use  of  certain  paraphrased 
words  for  the  various  quantities  encountered.  Table  I  shows 
these  paraphrased  words  also  its  equivalent  in  frequency  domain. 


TABLE  I 


Power  Cepstrum 
Domain 


Frequency 

Domain 


Definition 


Time 

Frequency 

Quefrency 

Cepstrum 

Analysis 

Cepstrum 


Complex 
nedonuil.it  ion 


Frequency 

Time 

Frequency 

Spectrum 

Analysis 


Spectrum 


Complex 
Demodula  t ion 


Procedure  of  summarizing, 
looking  at,  or  dissecting 
data 

A  dissection  of  the  vari¬ 
ance  of  time  series  into 
portions  associated  witli 
various  quef rencies 
(frequencies) 

A  shifting  of  frequency 
(frequency)  in  a  time 
series  by  multiplication 
by  sines,  and  by  cosines, 
of  a  carter  quefrequency 
(frequency),  followed  by 
smoothing  (and  sometimes 
decimation)  of  the  two 
res u 1 1 i ng  time  se r i es , 
which  can  he  regarded  as 
the  real  and  imaginary 
parts  of  a  complex  series 
(ana  logos  to  single-side¬ 
band  modulation) 
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TABLE  1  (Continued) 


Power  Cepstrum 
Domain 


Frequency 
Doma in 


Def inition 


Cross  Cepstrum 


Darius 
Gamnitude 
Lif tering 


Lifter 


Long  Pass  Lifter 


Lopar 


Cross  Spectrum 


Radius 

Magnitude 

Filtering 


F ilter 


A  dissection  of  the  com¬ 
mon  variation  of  two 
time  series,  into  portions 
associated  with  various 
quefrencies  (frequencies) 
which  takes  into  account, 
and  separates  both  syn¬ 
chronized  common  variance 
(covariance)  and  anti¬ 
synchronized  common  var¬ 
iance  (covariance)  and 
antisynchronized  common 
variation  (quadrature) 

Modulus  of  complex  number 

Modulus  of  complex  number 

A  transformation  of  one 
time  series  into  another 
which  a)  obeys  the 
superposition  rule  (is 
additive)  and  b)  is 
invariant  under  changes 
of  time  origin. 

A  device,  formula  or  pro¬ 
cess  for  making  such  a 
trans format ion 


High-Pass  Filter  One  which  passes  more 

rapidly  time  varving  com¬ 
ponents  more  readily 

Polar  Plot  or  coordinates  in 

terms  of  modulus  and 
ang  1  e 

The  number  of  a  t  ime 
ser i es  per  un it  time. 


Ouef rencv 


Frequency 


TABLE  1  (Continued 


Power  Ceps t rum 
Domain 


Rahmonic 


Repiod 


Saphe 


Frequency 

Domain 


Definition 


Harmonic 


Period 


Phase 


One  of  the  higher  que- 
frequency  (frequency)  com¬ 
ponents  generated  from  a 
sinusoidal  component  of 
time  series  by  a  non¬ 
linear  process  applied 
to  the  time  series  or  any 
equivalent  time  function 

The  amount  of  time 
required  for  one  cycle  of 
a  time  series 

Angular  displacement 
between  a  sinusoidal 
oscillation  and  a  ref¬ 
erence  cosinusoid  of  the 
same  quefrency  (frequency) 


Short-Pass  Lifter 


Los-Pass  Filter 


One  which  passes  rapidly 
varying  time  components 
more  readily 
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CEPSTRUM  DERIVATION  FOR  MULTIPLE  ECHOES 
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Extending  the  cepstrum  discussion  to  more  than  one  echo  or 
multi-path  signal  is  not  difficult.  We  limit  our  study  to  the 
case  of  an  original  signal  plus  two  echoes.  Suppose  z(t)  is  the 
signal  comprised  of  the  original  signal  y(t)  and  two  echoes,  then 
z(t)  can  be  written  as 

z(t)  =  y(t)  +  y  (t  -  t1)  +  a2  y  (t  -  t 2>  (B-l) 

taking  the  Fourier  transform  of  both  sides  in  (B-l) 

F(z (t) )  =  F(y(t)  +  a1  v  (t  -x^)  +  a2  y  (t  -t^)) 

=  F(y(t))  +  (XjF(y(t  -  T1))+a2F(y(t  -  x2>)  (B-2) 

where  F(  .  )  is  the  Fourier  transform. 

Calling  Fourier  transform  of  z(t)  and  y(t),  Z(cu)  and  Y(w) 
respectively  and  considering  time  shifting  property  of  the  Fourier 
transform  (B-2)  can  be  written  as 

Z(w)  =  Y (w)  +  u^y  (w)  e  ^  ‘  1  +  «7  y  (<»’)  e^’T  2  (B-3) 

or 

Z(w)  =  Y (w)  ( 1  +  e‘J“"l  +  ‘2e_1“n2)  (R_4) 


taking  absolute  value  of  both  sides 


.-.i'-'Si  i 2 


\7.  (<„)  |  =  I  Y(-.,)  (“j  (1  +  (]  e~’" 1  +  1 2  e  1  ''  2 ) 


(B-3) 


(B-3)  Is  comprised  of  two  parts 


'Y  (■■••) 


(B-h, 


and 


;  (l  + 


l  + 


1  9  e 


0  ' 


( B-7 ) 


j 
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expanding  ( B—  7 )  one  gets 


Id  +  «!  cos  <ox^  -  ja^  sin  an^  +  a ^  cos  an ^  -  jo^  sin  WT2^  I  ^ 

2 

=  Id  +  a1  cos  ut^  +  a2  cos  on  2)  -  j  (a^  sin  an  ^  +  02  sin  an 2)  | 

2  2 

=  1  +  +02  +  2  a^cos  idt^+  2a2  cos  01X2  +  a^oij  cos  o)(t2-Tj) 


+  0^2  cosu(t2  +  x^) 


(B-8) 


Lets  call 


2  2 

A  =  1  +  +  a2 


(B-9) 


then  (B-7)  can  be  written  as 


2a^  «2  ala2 

A  (1  +  cos  an^  +  2  —  cos  01X2+  — —  cosuCx^  -  x^) 


al°2 

+  — —  cos  01  ( x 2  +  i^)) 


now  taking  logarithm  of  both  sides  of  (B-5) 


(B-10) 


1  1 

2«n  |  Z  (w)  |  =  2  v"  n  |  Y  (<*>)  |  +  -  Vn  A  +  V.n  (  1  +  -7-  cos  on.  + 

Z  A  1 


"lu2  uirt2 

+  — —  cos  01  (  I  -i  +  — - —  cos  (.)  (  i  2+1  J  ) 


T  cos“’r2 


(B-ll) 


Note  that  and  .< ^  are  both  less  than  one  and  for  some  m, 1 ^  and 

12  we  can  have 


1  *■  *2  12  ']  '2 

1  •  — —  cos  1  ^  +  —7-  cos  M  l  2  +  -X —  cos  .,(  1  0-  :  ^  )  +  — :  cos  :  yt- 1  )  ;_  1 


(B-12) 


also  in(l  +  X)  can  be  expanded  as 


in  (  1  +  X)  =  X  -  ^  X2  +  ^  X3  -  t  X4 
l  3  4 


|  X  [  <  1  (B-13) 


Thus  up  to  the  first  approximation  (B-ll)  can  be  written  as 

1  2<X1  2<x2 

2  in  |  Z(w)  |  =  2 in  |  Y(w)  |  +  in  X  4 — —  cos  un  ^  4-  — —  cos  4- 

ala2 

cos  oj  +  — x —  cos  u  ^T2+Tp 


ala2 


(B-14) 


Taking  the  Fourier  transform  of  both  sides  of  (B-1A)  will  result 


Cp(z(t))  =  CEPSTRUM  =  F[2  in  1 2  (w)  j  =  F[2  in  |  Y (to)  |  ]  4- F[  in  X  ] 

2«i  2«2  rt1a2 

4-  F  [  — —  cosux  j  4-  — cos  tcT p  4- — —  cosw(t2-t^) 

“ia2 

4-  — —  cos  to  (t2+ti)  ( B- 1 5 ) 


or 

2^x 

C  ( z ( t ) )  =  CEPSTRUM  =  F[2  n|Y(w)i  )  4-  F[~  in  X  ]  4-  ■—  F[cos  oiT  ] 
p  2.  A  1 

2'*,  ,lllu2 

4-  —p  F  [cos  mi  ]  4-  — F[  cos  ui  (  r  i^)] 

(X 

+  ~~  Ffcos  ...(  i,-m ,)  ]  (B-16) 

(B-16)  has  several  terms,  the  first  part  will  give  a  regular  func¬ 
tion  not  necessarily  an  impulse,  the  second  term  gives  an  impulse 
at  the  origin  and  the  remaining  parts  would  give  impulses  at  the 
points  ip  i^,  '  9- '  j  and  1  a+  ]  on  t-*u‘  eepstrnl  plot  (Figure  C—  5  b). 

In  the  case  of  more  than  two  echoes,  the  procedure  is  the  same 
but  a  little  more  involved.  In  Appendix  C  we  sltow  eepst rum  plots 


for  the  cases  when  there  are  more  than  two  echoes. 


It  is  desired  to  calculate  the  cepstrum  for  a  signal  which 
is  a  combination  of  two  cosine  signals.  The  signals  have  started 
at  two  different  times  with  different  strength  and  although  both 
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doing  some  mathematical  manipulations  (C-3)  can  be  written  as 


F[f(t)]  =  — x - s—  [w(cos  A -cos  B  +  cos  C  -  cosD) 

(m  -  u>  Q) 

+  w  (cos  A-  cos  B  -  cos  C+  cos  D  ) 
o 

+  J 6o  (sin  A-  sin  B  +  sin  C  -  sin  D) 

+  Ju>  (sin  A-  sin  B  -  sin  C  +  sinD) 
o 

+  au>  (  cosR  -  cos  S  +  cos  T  -  cos  U  ) 

+  aco  (cosR  -  cos  S  -  cos  T  +  cos  U) 
o 

+  jau)(sinR  -  sin  S  +  sin  T  -  sinU) 

+  Jawo(sinR  -  sin  S  -  sin  T  +  sinU)] 

(C-4) 


where 


A 

= 

LOT 

+ 

o 

3 

1 

LO  T 
O 

B 

= 

COT 

o  -  w  T 

2  o  2 

+  LO 

C 

= 

0>T 

+  a)  t  - 

o 

LO  X 
O 

D 

= 

LOT 

«  +  T 

2  o  2 

-  fit 

R 

= 

WT 

-  a)  t  +  w  T 

o  o  3 


(0T  ■ 

-  O)  T, 

+  a)  t 

4 

o  4 

0 

COT 

-  CO  T  .  “  (0  T  ^ 

0  4 

o  J 

LOT  . 

+  U  I. 

-  (0  7 

4 

o  4 

0 

(C-5) 
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Equation  (C-4)  has  real  and  imaginary  parts,  to  find  its  magnitude 
both  real  and  imaginary  parts  should  be  squared.  Doing  the  pre¬ 
vious  operation  involves  lengthy  and  cumbersome calculat ions  that 
finally  would  result  the  following: 

2  2 

|  F(f  (t)  )  J  =  w  [4-2  cos  (A-B)  +  2  cos  (A-C)  -  2  cos  (A-D) 

-2  cos  (B-C)  -r  2  cos  (B-D)  -  2  cos(C-D)] 

2 

+  0>o  [4-2  cos  (A-B)  -  2  cos  (A-C)  +2  cos  (A-D) 

+  2  cos  (B-C)  -  2  cos  (B-D)  -  2cos(C-D)] 

2 

+  <u  [4-2  cos  (R-S)  -  2  cos(R-T)  +  2  cos  (R-U) 

+  2  cos  (S-T)  -  2  cos  (S-U)  -  2  cos  (T-U)  ] 

+  wwq[-4  cos  (A-B)  +  4  cos  (C-D)-4  cos  (R-S)  +4cos(T-U)] 

2 

+  nw  [  2  cos  (A-R)  -2  cos  (A-S)  +  2  cos  (A-T) 

-2  cos  (A-U)  -  2cos  (B-R)  +  2  cos(B-S) 

-2  cos  (B-T)  +  2cos(B-U)  +  2  cos(C-R) 

-2  cos  (C-S)  +  2  cos(C-T)  -  2  eos(C-U) 

-2  cos  (D-R)  +  2eos (D-S)  -2  cos(D-T) 


+2  cos  (D-D)  ] 
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2 

+  aoj^  [2  cos  (A-R)  -2  cos(A-S)  -  2  cos  (A-T) 

+  2cos (A-U)  -2  cos  (B-R)  +  2  cos(B-S) 

+  2cos (B-T)  -2  cos  (B-U)  -  2  cos(C-R) 

+  2  cos  (C-S)  +  2  cos  (C-T)  -  2  cos  (C-U) 

+  2  cos  (D-R)  -  2  cos  (D-S)  -  2  cos  (D-T) 

+  2  cos  (D-U)  ] 

+  oiwwo[4  cos  (A-R)  -  4cos(A-S)  -  4  cos  (B-R) 

+  4  cos(B-S)  -  2  cos  (C-T)  +  2  cos  (C-U) 

+  4  cos  (D-T)  -  4  cos(D-U) ]  (C-6) 

substituting  values  of  A,  B,  C  ...  in  each  term  of  (C-6)  will 
result  terms  that  are  completely  independent  of  io  and  also 
certain  terms  dependent  on  w  .  To  check  on  these,  lets  define 
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=  tu0('T2  +  T1  +  T  -  T3) 

* 1U  =  W0(~T2  +  T1  +  T4  “  V 

*11  =  0,o(-T2  +  -  T  +  T3) 

*12  =  Wo(~t2  +  T1  “  t4  +  t3> 

*13  =  Wo(“T  +  V 
*14  =  Wo(_2t  +  2t3} 


*15  =  Wo(-T  +  2t3  "  V 


*  =  w  (-2t,  +  2t„) 

16  °  A  3 


h  =  T  -  x  2 


Y2  =  T  -  t4 


Y3  r2  "t4 


( C—  7 ) 


after  substituting  (C-7)  into  (C-6)  then 


f  F  (  f  ( t )  )  f ' 


1  )  ^ 

— 2 - 2  2  f  1  +  Tf  "  2  cos  -(«)>!  +  *i ) 

(u)  -u  r 


-  2  cos  (uty i  +  i>3)  -  2  cos(-i.i>i  +  .}  ) 


2  cos  io  ( y  i  j)  -  2  cos(ioy2  + 


2  cos  (u)>2  +  *i5)  -  2  cos  (~uy g 


2  cos  ( toy ^  ~  ,'13)  -  2  cos  (u>  +  |fi) 
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UNCLASSIFIED  MSSU-EIRS-EE-Sl-J  AFOSR-TR-81-0603  NL 


-  2  cos  (wy2  +  4>g)  “  2  cos  (-wy1  +  i{ig) 


1  '  y9' 


+  2  cos  (wYj  +  4>^q)  -  2  cos  (-UY^  +  *11) 


+  2  cos  (wYg  +  4> 22 )  "  2  cos  (wy2  " 


-  2  cos  (wy2  -  <t> 3)  "  2  cos  ^-wVi  ~ 


+  2  cos  (wy3  “  “  2  cos  (-^Y-^  “  $9) 


+  2  cos  (wy3  ~  4>  10” 


+  ^  [C2  -  4  ui  QCOS (wy^  +  4>3)  +  4u>oCOs(uy3  _  4^) 


?  2 
4oi  woCOs(wy2  +  <j>33)  +  ^ o-  wq  cos  (uy2  ~  ^13^ 


-  4  a  u^cosCuYj  +  4>g)  ~  4  a  w  0cos  ]_  $9) 

+  4awocos(u>Y3  +  <J>jq)  +  2  a  w^cosCu^  ~  ^3) 
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with 

2  2 
C  =  ioQ  [4  -  2  cos  <{>2  -  2  cos  <|>^]  +  (aoi^)  [4-2  cos  -  2  cos 

2 

+  a  oi  o  [2  cos  4>^  —  2  cos  4>g]  =  CONSTANT 

2 

C^=  [4  +  2  cos  ^2  +  2  cos  <j>^]  +  a  [4  +  2  cos  <j>^  +  2  cos 
+  a  [  2  cos  (J>5  +  4  cos  4>?  +  2  cos  t)>9]  =  CONSTANT 

and 

C2  =  a  oio[4  cos  (|)5  -  2  cos  <j)7  ]  =  CONSTANT  (C-9) 

then  (C-8)  can  be  written  as 

I  F(f  <t) )  |  2  =  — 2~~2  ~2  11  +  G(ai)]  (C-10) 

(oi  -oi  ) 
o 

taking  logarithm  of  both  sides  of  (C-10)  one  gets 

2  , 

f.n | F ( f ( t ) )  |  =  in  — 2 - 2~2+  ?n  (1+  G(a)))  (C-ll) 

(o>  -01  ) 

o 

where  G(o>)  is  a  function  of  <o  . 

For  some  values  of  i->  ,  there  will  be  a  G(oi)  such  that 

|C(u»)  |  <_  1  and  then  ln(l  +  G(m))  can  be  expanded  as 

in(l  +  G  (o>) )  =  C(oi)  -  |(G(.o))2  +  ^  (G(oj))3  +  ...  (C-12) 

taking  the  first  term  of  this  expansion  (C-ll)  can  bo  written 

?n  |  F(f  (t)  )  |2  =  -2  '  n  (io2  -w  2)  +  G(u>)  (C- 13) 
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and  finally  cepstrum  is  found  when  the  Fourier  transform  of  both 
sides  of  (C-13)  is  calculated 

C  (f(t))  =  CEPSTRUM  =  F[Sln|F(f(t))  |2  =F[-2  in(w2-aio2)] +F[G(o)>] 

(C-14) 

Right-handside  of  (C-14)  is  comprised  of  two  parts,  first  part  which 
is  a  function  of  to  with  certain  Fourier  transform  and  the  second 
part,  Fourier  transform  of  G(w).  G(to)  is  a  set  of  cosine  wave¬ 
forms  such  as  cos  (toy  ^  +  <£  ^ ) ,  i  =  1,2,3,  J  =  1,2, ...  16  and  their 
Fourier  transforms  have  exponentially  weighted  impulses  at  points 
y^,Y2,  and  y^  .  Consequently  in  the  cepstrum  plot  peaks  will  be  at 
points  Y^  =  t— t 2 ,  y 2  =  t-t4»  and  ^3  =  T2~TU  (Figure  G-3b)  or  in 
the  other  word  there  is  a  peak  at  the  time  of  beginning  of  each 
signal  plus  a  peak  at  the  time  difference  of  two  starting  points 
for  the  signals.  Examaning  the  terms  in  (C-8)  will  reveal  that  the 
relative  magnitude  of  the  peak  at  the  point  Y3  =  T2~T4  is  small 
compared  to  those  for  starting  times.  Using  more  terms  of  the 
expansion  in  (C— 12)  would  result  in  a  term  which  causes  a  peak  at 
point  +  sunl  two  starting  points,  in  the  cepstrum  plots. 

The  computer  simulated  plots  for  this  case  of  our  study  can  be 
found  in  Figures  G-l  through  G-4. 


APPENDIX  D 

CEPSTRUM  DERIVATION  FOR  COSINE  TYPE  WAVEFORMS  WHEN 
THERE  IS  A  WINDOW  WITHIN  A  WINDOW 
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Figures  D-l  through  D-4  are  four  cosine  waveforms  which  have  dif¬ 
ferent  phases  and  strengths  .  The  equations  for  these  four  signals  are 

f^(t)  =  cos  0)ot  f 2 (t)  =  cos  WQ(t-T) 

f_(t)  =  a  cos  u)  (t-v)  f(t)  =  cos  u  (t-i)+acosw  (t-v)  (D-l) 

3  °  o  o 

It  is  obvious  from  (D-l)  that  f(t)  is  the  combination  of  f ^ ( t)  and 

f^(t).  These  four  signals  have  all  started  at  times  t^  234 

seconds  with  t^  equal  to  t^  or  t^  depending  which  is  the  smaller  one. 

Now  suppose  we  want  to  study  the  cepstrum  for  f(t)  while  it 
is  observed  through  a  time  window  having  length  X-y  seconds  which 
is  defined  by 

u(t-Y)  -  u(t-X)  (D-2) 

as  is  shown  in  Figure  D-4  with  dotted  lines.  The  procedure  to  find 
the  cepstrum  is  finding  the  Fourier  transform  of  the  signal,  the 
squared  spectrum,  the  log-squared  spectrum,  the  Fourier  transform, 
and  finally  the  spectrum  again.  Now  with  regard  to  all  above 
operations  we  proceed  to  do  the  calculations. 

F[(f(t))]  =  /  [cosm  (t-r)  +  a  cos  iii  (t-v)] 

o  o 

[u(t-y)  -u(t-X)]e''ultdt  (D-3) 

Term  [u(t-v)  -  u(t-X)]  is  a  combination  of  two  unit  step  functions 
which  cause  uhe  limits  of  the  integration  to  be  y  and  X.  Mow  with 
these  new  limits  (P-3)  becomes 
.X 

Ff  f  (  t )  1  =  j  f  COS  'n  (t-i) 

Y  0 

+  1  cost  (t  -  )  ]  e  ^  tdt 


(D-4) 
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using  the  Euler's  identity  for  cosine  (D-4)  can  be  written  as 


F[f(t)]  =  / 


x  r 


i  [eja,o(t-T)  +  e-j“o(t-T)] 


+  \  [ejwo(t"v)  +  e-ja)o(t"v)]]  e_J“tdt 


(D-5) 


—•jut 

multiplying  e  into  the  bracket  and  taking  all  variables  not 
dependent  on  t  out  of  integral  one  gets 


_  r  /  \  i  1  -iw  T  r ^  -j  (lx)-UJ  )t  , 

F[f  (t)  ]  =  ^  e  J  o  J  e  J  o'  dt 

Y 


,1  iw-r  r  X  -  i  (oo+ai  )  t 
+  -=■  e  o  J  e  o 


dt 


,  1  -jai  v  cA  -j  (uj-di  )t  . 

+  o  Je  o  dt 

Y 

,1  iu  V  r  ^  -j  (w+u)  )  t  . 

+  j  aeJ  o  J  e  J  o  dt 


(D-6) 


(D-6)  is  comprised  of  four  parts  each  having  a  term  such  as 


r  p.t  1  ,  ax  &y 

j  e  dt  =  —  (  e-e) 


(D-7) 


When  these  four  terms  are  integrated  and  combined  the  result  is 

...  1  r  -  ]  (<ii  t+ojA-ojA)  -  j  (<u  i+wy-di  i) 

F[  f  (t)  ]  =-T^r - 1  fe  0  °  "e  o  ° 


j  2  (u)— io  ) 


,  -  i  ( u)  y  +  in)  -di  A  )  -j  (m  o  +  ioy  -  hi  Y  )  ] 

+  a  e  J  o  o  ->e  o  o 


1  .  -if-io  t  +  i.i  '  +  (<!  '■)  -  j  (-hi  i  +  ..)  +  (i  y) 

-  1  e  o  o  -e  o  o 


i  2  (io+h i  ) 

o 


—  j  (— o>  '  +iii\  +  io  ')  ~i(~:o  +  h'Y  +  io  y)] 

+  y  e  J  o  o  - 1  e  o  o 


(I)-8) 
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with  (D-8)  having  real  and  imaginary  parts. 

The  spectrum  of  a  signal  is  found  by  squaring  the  real  and 
imaginary  parts  and  then  talcing  its  square  root.  In  doing  this 
after  some  mathematical  manipulations  one  gets 

2  1  1 

I  F(f  (t)  )  |  =  [— 2 - 1  [  1  +  -a  (COS  to  (t-  y)  +  3  ot  cos  to  ( ■>  -y  ) 

to  -10  *■  8  0  0 

o 

2  2 

-  (1+a)  cos  (to-to  )<J>  -  (1+a  )  cos  (to+to  )cj>  -  acos(to  (t-y)+(io-io  )d>) 

o  o  o  o 

-  a  cos(to  (x-y)  -  (to-to  )<f>)  -  a  cos(-to  (t-y)  +  (to+to  )d>) 

O  O  0  0 

-a  cos  (to  (t-y)  +  (to+to  )40]  (D-9) 

o  o 

where 

<t>  =  A  -  y 

To  find  the  cepstrum  logarithm  of  (D-9)  should  be  taken. 

As  is  shown  previously  £n(l  +  X)  when  ]x|  <_  1  can  be 
/•expanded  to 

en(l  +  X)  =  X  -  |  X2  +  j  X3-  j  XA  +  ...  (D-10) 

also 

Vn  AB  =  in  A  +  in  B  (D-ll) 

Knowing  t  <_  1  and  also  |cos  x|  •_  1  ,  it  will  be  justified 

to  have  assumed  that  the  term  in  the  parenthesis  in  (D-9)  is  less 
than  one  and  with  this  (D-9)  can  be  written  as 
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in  |  F(f  (t))  j  ^  =  In  {— — ^ — r  ]  +  i(c os  u  (  r-y)  +  3a  cosu>  (t  -y) 

z  z  o  o  o 

a)  *"Co 

o 

2  2 

-  (1  +  a  )  cos  (w-u>o)iJ>  -  (1  +  a  )cos(w+wo)  $  -  a  cos(u>o(t-y)+(u>~u)o)<}>) 

-  a  cos  (wQ(x-y)  -  (m-mo)4>)  -  a  cos  (-w^ft-y)  +  (co+w^)  <p  ) 

-  a  cos  (to  (t-y)  +  <oj-Koo) 4> )  (D-12) 

The  final  step  in  finding  cepstrum  of  f(t)  function  is  to  find 
the  Fourier  transform  of  (D-12)  again  with  this  in  mind  that  the 
variable  of  integration,  representing  time,  here  is  u  instead  of 
t.  With  a  glance  at  (D-12)  presence  of  several  cosines  such  as 
coswo(t-y)  =  CONSTANT  and  cos(A  +  w<j>)  with  A  as  a  constant 
phase  is  obvious.  The  Fourier  transforms  of  these  terms  is 

B  =  A  CONSTANT-* - *  B5(f) 

cos  (w  t+0)+ — >•  4[e^  (f-f  r)+e  ^  d(f+f  )  ] 
c  c  c 

(D-13) 

which  shows  impulses  at  t  =  0  second  and  t  =  f  seconds. 

c 

Translating  terms  such  as  these  into  the  cepstrum  plot  would  give 

peaks  at  t  =  0  and  t  -  \  seconds  which  proves  the  existence 

of  peaks  in  the  cepstrum  plot  only  at  points  K ( A  —  y ) ,  ,  _ 

K-i  . 

(length  of  time  window)  and  its  harmonics. 


APPENDIX  E 


CEPSTRUM  DERIVATION  FOR  DAMPED  EXPONENTIAL  TYPE  SIGNALS 
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Figure  (E-l)  shows  three  signals  which  are  exponentially 
damped.  Figure  (E-lc)  is  the  combination  of  Figures  (E-la)  and 
(E-lb) .  The  equations  for  these  signals  are 

fa(t)  =  (t-Tx)  e“C(t_Tl)  u  (t-Tx) 

f2(t)  =A(t-x2)  e~C(t_V  u  (t-T2) 

f(t)  =  fx(t)  +  f2(t) 

=  (t-ij)  e  c^t  Tl^u(t-r1)  + A  (t-r2)  e  c^t  T2^u(t-T2) 

(E-l) 

Fourier  transforming  f(t)  when  it  is  observed  from  time  t=0 
to  time  t=r  seconds  would  result 

F  [  f  ( t )  ]  = 

+ 

(E-2)  can  be  b 


/  f(t)  e  ^W<dt 


/  (t-T1)^C(t'r])e-j'Jldt 


;  ,  ,  .  -C(t-T„)  , 

J  A ( t- c  )  e  2  e  J  dt 


(E-2) 


roken  in  four  parts  as  follows 
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F[  f  (t)  ]  =  /T  te-(C+ju))teCTldt 
T1 

-  Ti/Te-(C+jco)teCll  dt 
T1 

+  kj  \  e'(C+ju)teCl2dt 
T2 


-At 


2 


f  e“(C+jlj)teCT2dt 

t2 


(E-3) 


Lets  have  these  new  notations  for  the  simplicity  reason 


(1)  =  /Tte'(C+ju)teCTldt 

T1 

(2)  =  -x1  /  e-(C+jw)teCTldt 

T1 

(3)  =  A  /T  te-(C+'1w)TeCT2dt 

T  2 

T 

(4)  =  -At 2  /-e  (C+ju>) t  eCT2(lt 

r2 

So  (E-3)  can  be  written  as 


Fff(t)]  =  (1)  +  (2)  +  (3)  +  (4) 


(E-5) 


To  find  the  Fourier  transform,  each  term  is  calculated  separately  and 
then  the  results  will  be  added  to  give  the  Fourier  transform  of  f(t) 


—  1— ■■ 


V 
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and  finally 


(4)  = 


AKC 

-  )  AKu 

COS  OJT  + 

-KAco 

-  jAKC 

7 

,  2 

2 

2 

C 

{ 

4*  a) 

C 

+  10 

sm  cot 


-At2C  +  j  Acot  2 

cos  cot2  + 

At2co  +  JAt2C 

J2  ,  2 

2  2 

C  +  to 

C  +  to 

.  J 

sm  cot. 


(E-12) 


with 


K  =  T2eC(T2"T) 


(E-13) 


The  Fourier  transform  of  signal  f(t)  is  sum  of  (1),  (2),  (3),  and 
(4).  Doing  summation  operation  one  gets 


F[  f  (t)  ]  =  (D1  +  jD2)  cos  cot  +  (D2  -  ;jDj)  sin  cot 

+  (D^  +  ,TD4)  cos  cot  (0^  -  jD^)  sin  cot^ 

+  (D,-  +  jD^)  COSWT2+  (D^  -  jD^)  sin  cot2 

where 

=  Mo  -  NfcSco  +  Pu  -  0  Bio  +  R^  +  R2 

D2  =  -6M10  -  Nu  -  PB10  -  Qn  -  T  -  T2 

D ^  =  -  |  M  +  N  A  <0  -  SI 

D.  =  M  +  NY  +  U1 
4 

D5  =  -Pv  +  Oou  -  S2 


(E-14) 


D,  =  P  u>  +  vO  +  1-2 
o 


(E-15) 
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and 


R1  " 


KC 


„2.  2 

C  +o) 


AKC 

2  2 

C  +w 


S1  " 


T1C 

c2+.2 


S2  = 


At2C 


2  2 
C  +  w 


Ti  = 


Hw 


2  '  2 
C  +  to 


T  = 
2 


AKw 

2  2 
C  +  u) 


~2  2 
C  +w 


U2  " 


At  2^ 

T2~  2 

C  +o) 


2  2 


M  =  eCTl  -4^ 


2  2  2 

(C  +co  ) 


N  =  e 


Cri  2ojC 


2  2  2 

(C  +<o  ) 


2  2 

„  .  C  T  9  C  -0.1 

P  =  A  e  2  — - — y 

C  +,o“ 


Q  =  A^[2 


„2  2 

C  -Ho 


(E- 


Next  stop  in  calculating  the  cepstrum  is  of  finding  square  of 
absolute  value  of  (E-3).  After  some  mathematical  manipulations 
the  final  result  is 


j  F  ( f  ( t )  ) 


9  p  a  a  9  “>  9 

=  (D1  +  D2  +  D3  +  D4  +  D5“  +  Dfe  ) 

+  (  2D  I)„  +  2i)„D.  )  cos  io ( i  —  i  ,  ) 

1  J  Z  4  I 

+  (2f>2n3  -  2n  n4)  sinw(r-i  j) 

+  (2D,  D.  +  2D_  0r  )  cos  <4i(r  -r  ) 

4  6  3  3  2  1 


Li 


16) 


+  (20,0,.  +  20„0  )  cos  .i(  i-  O 
15  2  b  2 
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+  (2D^D^  -  2D2D^)  sin  us(t  2~t ) 


+  (2D^D5  -  2D^D^)  sincj(i1-i2) 


(E-17) 


To  find  the  cepstrum  logarithm  of  (E-17)  should  be  taken. 

As  in  the  previous  case  (E-17)  can  be  written  as  1  +  G(w)  and 
for  some  u  ,  Hn  (1  +  G(w)  =  G(to).  Then  when  the  Fourier  transform 
is  taken,  there  should  be  weighted  impulses  at  ,  t ^ ,  t 2~ri>  ^ut 
as  we  will  show  with  an  example  the  strength  of  these  impulses 
depends  strictly  on  the  value  of  the  damping  coefficient  C. 


Example:  Let 


C  =  30.0 


t  =  . 5  sec 


t  =  .07  sec 
i2  =  .12  sec 
A  =  .  6 


2  2  2 

D  =  -1.47xlO~4  +  2 . 77x1 0~ 4  — +  9. 8lxl0~7  — -x 

C  +u>  (C  +ui  )  C  +i.i 


D  =  -4 . 61xl0_6  }  +  8.45x10  3  y 

L  cz+..r  c^+i.r 


2  2  2 

D  =  3  i  +  2  j 

J  1  9  9?  ^  *  -1  9  ? 

(C"+</)  c  +./ 


D  =  0.07  -(4 -186  ~~  +  .07 

(€*•+</)  c  +.,.  c  +..r 
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1. 56u)2  +  2.76  C2  -  8.64 

2  2 
C  +u 


2  2 

D  =  7 . 2xl0-2  -166  -2-5- 

6  C2+w2  C2+w2 


f  Substituting  values 

\ 


through  D, 
o 


into  (E— 17)  will  show  that  the 


only  term  with  noticeable  strength  is  the  coefficient  of 


cos  to  ( J2~ri )  or  the  only  visible  peak  in  the  cepstrum  plot 
is  the  peak  at  the  time  difference  t2_ti  *  For  colnPuter 
simulated  results  see  Figures  G-ll  through  G-14. 


* 
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A  SYNOPSIS  AND  LISTING  OF  DIGITAL  COMPUTER  PROGRAMS 


i 
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A  SYNOPISIS  OF  DIGITAL  COMPUTER  PROGRAMS  AND  THEIR  LISTING 

The  Univac  1108  computer  at  Mississippi  State  University  Com¬ 
puting  Center  was  used  for  the  development  of  all  runs  and  simula¬ 
tion  plots  in  this  study.  The  programs  are  all  written  in  FORTRAN 
.  IV. 

The  programs  are  all  written  in  a  machine  independent  form, 
such  that  they  can  be  run  in  any  type  of  computer  with  very  little 
change.  Effort  has  been  done  to  write  the  programs  as  clear  and 
as  sufficient  as  possible  so  the  computer  run  time  and  consequently 
the  cost  would  be  minimized.  There  are  enough  comments  in  each 
program  to  clarify  any  operation  in  the  programs.  At  this  point 
it  is  appropriate  to  give  a  brief  summary  of  each  program. 

MAIN  PROGRAM 

The  main  program  reads  and  writes  the  initial  data  input  also 
calculates  the  power  spectrum  after  taking  the  Fourier  transform 
for  the  first  time,  then  it  finds  the  logarithm  of  the  square  of 
the  power  spectrum  to  send  it  to  the  subroutine  FFT  for  Fourier 
transforming  again,  it  then  calculates  the  Cepstrum. 

The  number  of  sample  points  for  the  original  signal  is  2048 
samples/second  and  the  total  period  of  sampling  is  0.5  seconds.  We 
tapered  the  first  and  last  20  points  at  each  end  of  the  cepstrum 
data.  This  tapering  operation  involves  multiplication  of  the  i*-*1 
value  from  each  end  by  1/2  [1  -cos  (i/20)],  thus  making  a  smooth 
transition  from  values  multiplied  by  zero  to  those  multiplied  by 
unity.  Finally  we  plot  the  cepstrum  versus  time  to  observe  its  peaks. 
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SUBROUTINE  SIGNAL (Y)  (Cosine) 

This  subroutine  generates  few  cosine  waveforms  with  different 
strength  and  phase  but  same  frequency.  The  main  signal  is 

y  (t)  =  cos  (  2 irf  (t— t^)  u  (t-t^) 

and  the  echoes  are 

y ( t)  =  ol^  cos  (2  'i  f(t-T^))  u  (t-iO 

a .  <  1  . 
r  — 

where  f  is  the  frequency  and  is  the  epoch  or  echo  delay  time 
the  composite  signal  is  the  one  whose  cepstrum  is  sought. 

SUBROUTINE  SIGNAL  (Y)  (Cosine  plus  Noise) 

This  subroutine  is  like  the  previous  one  except  that  it  mixes 
the  white  Gaussian  noise  with  the  composite  signal  to  make  a  noisy 
signal.  The  complex  array  Y  is  generated  and  sent  to  the  main 
program.  Statistical  property  of  the  noise  is  investigated, 
signal  to  noise  ratio  for  the  signal  is  also  calculated. 

SUBROUTINE  SIGNAL  (Y)  (Damped  Cosine) 

Subroutine  signal  for  damped  cosine  will  generate  a  f ew 
damped  cosine  waveforms  of  the  type 

v  i  ( t )  =  i .  cos  (2"f(t-T^))e  ^  ^ t  i^u(t-i.) 

a.  =  •  1  . 

1 


each  with  different  strength.  Parameters  f  and  x^  are  as  in  the 
previous  case  and  C  is  the  damping  coefficient.  All  variables  can 
get  different  values  furnished  through  the  main  program. 

SUBROUTINE  SIGNAL  (Y)  (Damped  Exponential) 

In  this  subroutine  few  damping  exponential  signals  are 
generated  each  starting  at  different  times  and  having  different 
strength.  The  signals  are 

Y^Ct)  =  <*^(t-x^)  e  Ti^u  (t-x^) 


Again  t.  s  are  the  epoch  or  echo  delay  times.  These  signals  will 
be  added  together  to  give  the  composite  signal. 

SUBROUTINE  WINDOW  (I,N,W) 

This  subroutine  does  the  windowing  process  based  on  Kaiser- 
Bessel  equations  with  ft  =  8.0.  I  is  the  sample  point,  N  is  the 
total  number  of  samples  and  W  is  the  weight  which  is  calculated 
in  this  subroutine. 

SUBROUTINE  EFT (X,M) 

This  subroutine  will  calculate  the  Fourier  transform  of  a 

M 

complex  array  X  with  total  number  of  data  to  be  2  .  The  routine  is 
a  Fast  Fourier  Transformation  based  on  Cooly-I'ukey  algorithm. 
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The  complex  array  X  is  returned  as  the  Fourier  transform  of  the 

input  array  X.  The  maximum  value  of  M  is  10  which  can  be  increased  | 

if  the  respected  dimensions  in  all  routines  of  this  report  are 
changed . 

SUBROUTINE  SMOOTH  (N,  Y)  | 

Smooth  subroutine  does  smoothing  process  by  Hanning  method, 

i 

N  is  the  total  number  of  points  and  Y  is  a  complex  array  of  data  to  > 

be  smoothed.  After  the  process  is  done  array  Y,  the  smoothed  one, 
would  return. 

SUBROUTINE  RANSET  (MAXINT,  NSTRT) 

Subroutine  RANSET  in  coorperation  with  subroutine  URAND  con¬ 
stitute  a  machine  independent  random  number  generator  having 
specially  good  statistical  and  correlation  properties.  MAXINT  is 
a  very  big  integer  depending  on  size  of  the  computer  (for  example 
32000  for  a  PDP-11  computer)  and  NSTRT  a  small  integer  usually  0. 

REAL  FUNCTION  URAND  (FRAN) 

This  real  function  subprogram  has  been  used  in  cooperation 
with  subroutine  RANSET  to  generate  the  white  Gaussian  noise.  FRAN 
is  a  dummy  variable  which  can  be  anything.  The  necessary  parameters 
needed  in  this  program  are  furnished  through  a  common  statement  from 
RANSET  subroutine. 
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SUBROUTINE  SHEKL  (B , A , NP , XNAMF. , YNAME , EMI t  EX1 , IC 1 ) 

Subroutine  SHEKL  is  a  plotting  subroutine  which  plots  array  A 
versus  array  B.  It  is  a  machine  independet  plotting  routine  with 
NP  being  the  total  number  of  points  to  be  plotted  and  not  necessarily 
the  total  number  of  data  in  arrays  A  and  B.  XNAME  and  YNAME  are 
two  16  character  axis  title  for  X  and  Y  axis.  EMI  and  EX1  are  the 
minimum  and  the  maximum  given  to  the  Y-axis  if  desired  to  do  so  and 
finally  IC1  is  set  equal  to  zero  if  no  desire  for  setting  maximum 
and  minimum  for  Y  axis  and  is  set  equal  to  one  otherwise. 

SUBROUTINE  RPRD 


This  is  the  only  machine  dependent  routine  used  in  the  entire 
study.  It  is  used  to  plot  the  Cepstrum  and  is  fed  for  its  variables 
through  common  statements  by  the  main  program  and  by  the  signal 
subroutine.  SIGCMP  is  a  common  statement  bring  variables  from 
the  signal  subroutine  and  CMPLOO,  CMPLOl  and  CMPL02  are  three 
common  statements  responsible  for  transfer  of  variables  from  the 
main  program  to  the  subroutine. 


Ill 


A  LISTING  OF  FORTRAN  PROGRAMS 

The  following  are  programs  and  subroutines  explained  completely 
in  previous  pages. 
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c 

c 

1 

2 

C 

C 

C 


*  *  *  * 
*  *  *  * 


#  *  *  * 


COMPLEX  Y  (10  2  0 

DIMENS  ION  AMAGC1024),XNAME  C4)  ,YNAME(4  >,PC102O 
COMMON/CMPL2 /NP,T  ,F  ,ALf 0  ,  ALF1  .COEf  ,TAU(3> 

COMMON  /  0  AT  AA  / HN AME  (A)  ,  V  NA  f*E  (4  )  ,Q(1u24) 

DATA  XNAME/4hFREQ,<,hUEN-,4hCY  ,4H  / 

DATA  YNAME/4HMAGN,4HITU-,4HDE  ,4H  / 

DATA  HNAME/4HTIME  ,4h  S£C,4h0ND  ,4H  / 

DATA  VNAME/4HSI6N.4HAL  ,  4H  ,4H  / 

IWNDW  IS  f OR  WINDOWING  IN  TIME  IWNDW  1  YES  0  NO 
IWNDF  IS  FOR  WINDOWING  IN  FREQUENCY  1  YES  0  NO 
FORMAT ( 1 H 1 ) 

FORMAT (3 (/>> 


SET  ALL  CONSTANTS  USED  IN  THE  PROGRAM 


IC  =  2 
£M=  1 , 

£  X  =  1  . 

JJ  =  1 

PI=3. 141592654 
PI20=PI/20. 

C 

C 

C 

C****  READ  DATA  FROM  DATA  CARD 
C 

REA0C5  ,5 )  M,ND,IWNDW,IUNDF,T,F,ALFO,ALF1,COEF 
5  FORMAT  ^4 12  ,5  *10 .5  ) 

REA0C5.10)  (TAU (1 )  ,  1  =  1 ,ND  ) 

10  FORMAT  (8  F 1 0, 5) 

C 


c 

c 

C****  PRINT  THE  READ  DATA 
C 

WRITE  (6,1  ) 

WRITE (6,12)  M,T,NO ,IWNOW,r«NOf ,F, A LF0.ALF1, COE F 
12  FORMAT (23x,2hm=, I2/23X ,2HT= ,F10.5/21X ,4H  ND=,I2/ 

«19x,6H1wnDW=,I1/19x,6H1wN0F=,I1/23x,2hf=,F10.5/ 
*20X,5HALFO=,F10.5/20X,5HALF1=,F10.5/20X,5HCOEF=, 

*  f  1  0  «  5  ) 

DO  20  1  =  1,  ND 

WRITE  (6.15  >  I.TAUCI) 

15  FORMAT (17X,4hTAU(,12»2H)=»F10»5) 

20  CONTINUE 

C 

c 

NP  =  2*  *M 

IFC JJ  .  N  E  •  1  )  GO  TO  26 

C 

C  *  *  *  *  CALCULATE  FREQUENCY  FOR  THE  FOURIER  TRANSFORM 
c***«  NOTE  that  WE  ARE  USING  FFT  AND  THE  FREQUENCY  IS 
C  *  *  *  *  (NUUFSJ/NP  WHERE  FS=1/TS  AND  T  S  =  S  AMP  L  2  NG  FREQUENCY 
TS=T/NP 
f  $  =  1  /  T  S 


FSS=FS/NP 


TN=T/NP 

P(1)=0. 

NP1  =NP-1 
DO  25  1=1, NP1 
Q(I)=TN*I 
NU= 1*1 
P(NU)-I *FSS 
25  CONTINUE 
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C 

C 

Q(NP)=NP*TN 

C 

C 

c 

O***  GENERATION  OF  THE  SIGNAL 
CALL  SIGNAL(Y) 

C 

c 

26  CONTINUE 

IF(  IWNDW.EQ.O)  GO  TO  40 
C 

C‘***  DO  THE  WINDOWING  PROCESS 
C 

DO  30  1=1, NP 

CALL  WINDOWCI ,NP,W) 

X  =  W*REAL  (Y  (I  )  > 

Y(I)=CMPLX(X,AIMAG(Y(I))) 

30  CONTINUE 

40  CONTINUE 

C 

C  *  **  *  FIND  THE  FIRST  FOURIER  TRANSFORM 

C 

CALL  FFTCY.M) 

C 

C 

C 

C****  FIND  THE  LOG  MAGNITUDE  ALSO  MAKE  PREPARATION  FOR 
c***<  TAKING  THE  SECOND  FOURIER  TRANSFORM 
C 

DO  60  1  =  1, NP 

AMAG(I)=REAL(YCI))**2'*AIMAG(VCI))**2 
IFCAMAG(I).EQ.O.O)  GO  TO  50 
X=ALOG10(AMAG(I)> 

V(I)  =  CMPLX  (X  , 0.0) 

GO  TO  60 
50  CONTINUE 

X=-1 .0  El  0 
Y(1  )  =  CMPLX  CX ,0.0) 

60  CONTINUE 

c 

c 

c 

C  *  ** •  PLOT  THE  POWER  SPECTRUM 
C 

WRITEC6.2) 

WRITEI6.65) 

65  FORMAT (lOX  ,34HTH1S  IS  GRAPH  OF  THE  FIRST  FOURIER, 

*33H  TRANSFORM  OF  THE  COMPOSIT  SIGNAL) 

WRITE  (6,2) 

C 

NP  =  NP / 2 

CALL  SHEKL(Q,AMAG,NP,XNAME,YNAME,EM,EX,IC) 

NP=2*NP 

C 

IF(IWNDF.EQ.O)  GO  TO  70 
C 

C*  ***  DO  THE  SMOOTHING 
C 

CALL  SMOOT H(NP , Y) 

C 

C 

70  CONTINUE 
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114 


C 

C  *  **  *  FIND  THE  FOURIER  TRANSFORM  OF  THE  LOG  SPECTRUM 

CALL  F  FT  <  Y  ,M  ) 

C 

C 

c 

C*“*  CALCULATE  THE  POWER  CEPSTRUM 
C 

WRITE(6,1  ) 

NP  =  N P/2 
DO  90  1=1, NP 

AM  A  G  (  I  ) =S  Q  R  T  (  R  E  A  L  (  V  (  I  )  )  *  *  2*  AI  MAG  £  Y  (  I  )  )  *  *2  ) 
1FC1.GT.20)  GO  TO  90 
WRITE  (6,80)  Q(I>,AMAG(I) 

80  FORMAT  (2 C10X  ,620.6)) 

90  CONTINUE 

C 

C 

c 

C  * » *  *  DO  SMOOTHING  OF  THE  £N0  POINTS  PY  MULTIPLYING 
C*  ***  THE  FIRST  20  POINTS  BY  .5  * ( 1 .-COS < 2 . *PI * 1/ 2 0. > ) 
C*“*  1  =  1,2, . 20 


WRITE  f 6 , 1  ) 

WRITE (6.120) 

120  FORMAT (10X.25HTHIS  IS  THE  CEPSTRUM  PLOT) 
WRITE(6,2) 


CALL  SHEKL(Q,AMAG,NP,XNAME,YNAME,EM,EX,1C) 


STOP 

END 


MAIN  PROGRAM  (Continued) 
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*  *  *  * 
*  •  *  * 
*  *  *  * 
•  *  *  * 
*  *  *  * 
*  *  *  * 
*  *  *  * 


*  *  *  * 


4*0* 


0 


4  4  4  4 


SUBROUTINE  SIGNAL(Y) 

COMPLEX  V  < 1 0  2  4 ) 

DIMENSION  AMAGC1024),AMGAC1024>,CMAGC1024) 
DIMENSION  DMAGC1024) 

COMMON/ CMP L2/NP,T  ,F,ALF0,ALF1,C0EF,TAUC3> 

COMMON /D  AT  AA /HNAME  C4),VNAM£C4),QC1024) 

THIS  SUBROUTINE  WILL  GENERATE  FEW  COSINE  WAVEFORMS 
WITH  THE  SAME  FREQUENCY  BUT  DIFFERENT  STRENGTH 
AND  ARRIVING  AT  DIFFERENT  TIMES 

TAU ( I )  1*1,2,..  IS  THE  DELAY  IN  SECONDS 
F  is  the  frequency  in  herti 


FORMAT  C 1 H 1 ) 

FORMAT (3 ( / ) ) 

SET  THE  CONSTANTS 


EM  =  1  . 

E  X  =  1  . 

I C  =  2 

PI=3. 1 41 592654 

PI2=2.*PI 

P12  F  =  P 12  *  f 

TN=T/NP 

SS*0.0 

T1  =  .2 

T2  = .35 

CS3=0. 


SIGNAL  GENERATION 

DO  20  1=1, NP 
TM  =  TN*  1 
TM1  =TM-TAU  Cl  ) 
TM2=TM-TAUC2) 

TM3=TM-T  AU  (3  ) 

TM1 P  =  P 12  F*TM  T 
TM2P=P12F*TM2 
TM3P=PI2F*T"3 

CS1=C0S(TM1P) 

CS2  =AL  FO*  COS  (TM2P) 
CS3=ALF1*C0S  (TM3P) 

IFCTM .LT.TAU  Cl ))CS1 =0.0 

IFCTM.LT .TAU C?)  )  CS2=0.0 

IFC TM.LT .TAU C 3)  )  CS3  =  0. 

CS=CS1 *CS24CS3 

AMAGC 1 )  =  CS  1 

AMG AC  I ) =C  S  2 

DMAGC  J  )  =C  S  3 

CMAGC I >=CS 

YC I)=CMPLX  CCS  ,SS> 

CONTINUE 


PLOT  THE  SIGNALS 

MP*  300 
WR1TEC6.1  ) 
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WRITE (6,30) 

30  FORMAT ( 1  0  X  ,33HTH1S  IS  G RA PH  OF  THE  FIRST  SIGNAL) 
WRITE(6,2) 

CALL  $HEKL(Q,AMAG,MP,XNAM£#YNAME,EM,EX,IC) 

WRITE  (6,1  ) 

WRITE  (6.40) 

40  FORMAT (10X,34HTHIS  IS  GRAPH  OF  THE  SECOND  SIGNAL) 

WR1 TE  (6,2) 

CALL  SHEKL (Q,AMGA,MP,XNAME,YNAME, EM,  EXtIC) 
IFXTAU(3).EQ.0.)  GO  TO  60 
WRITE(6,1  ) 

WRITE (6.50) 

50  FORMAT (10X  ,33HTHIS  IS  GRAPH  OF  THE  THIRD  SIGNAL) 
WRITE(6,2) 

call  SHEKL(Q,DMAG,MP,XNAKE,VNAME,EM,EX,1C) 

60  CONTINUE 

WRI TE (6 ,1  ) 

WR1TE(6,70) 

70  FORMAT  (10X  .36HTHIS  IS  GRAPH  OF  THE  COMPOSIT  SIGNAL) 
WRIT£(6,2) 

CALL  SHEKL(Q , CM AG , MP , XN AM E , YN AME , E M , E X , I C ) 

WRITEC6.1  ) 

RETURN 

END 


SUBROUTINE  SIGNAL  (Cosine)  (Continued) 
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C  *  *  *  • 
C  •  •  •  • 

c  *  *  *  * 
c**«* 
*  .  »  . 
*  *  *  * 


SUBROUTINE  S IGNAL(T) 

COMPLEX  r<T024) .ANC1024) 

DIMENSION  AMA6 11024 > ,BM AG  £1  02  4  >  ,  C  M  AG  <  1  024  >  ,  AN  S  1 1024  ) 
DIMENSION  DM  A  G  £  1  02  4  ) 

COMMON/CM2/NP,TfF,ALFO,ALF1.COEF.TAUC3),FSNR,M 
COMMON  /  DAT  A  A  /HNAME  £  4  )  ,  V  N  A  ME  (4  )  ,  0  (  1  024  ) 


THIS  SUBROUTINE  HILL  GENERATE  FEW  COSINE  WAVEFORM 
WITH  THE  SAME  FREQUENCY  BUT  DIFFERENT  STRENGTH 
AND  ARRIIV1NG  AT  DIFFERENT  TIMEStAlSO  IT 
GENERATES  AND  MIXES  NOISE  WITH  THE  SIGNALS 

FORMAT (1  Hi ) 

FORMAT  £3  £/>) 


****  SET  ALL  CONSTANTS  USED  IN  THE  PROGRAM 


PI=3. 141592654 
PI2  =2 • *P  I 
PI2  F -P 12 * F 
TN-T/NP 
SS=0.0 
SSN  =0.0 
E  M  =  1 . 

E  X*  1  • 

I C  =  2 
SUM =0  .0 
SUM2  =  0 .0 
CS3 -0 . 0 


***  INITIATE  THE  NOISE  GENERATOR 
CALL  RANSET(1E*10,0) 


****  GENERATION  OF  THE  SIGNAL 

DO  20  1*1  ,NP 

TM  =  TN  *  I 

TM1 =TM-T AU  £1  ) 

TM2  =TM-T  AU £2  ) 

TM  3  =T  M  -T  A  U  £3  ) 
TM1P*PI2F*TM1 
TM2  P  =  PI2F*TM2 
TM3P  =  P 12  F  *T«3 
CS1 *C0S  £ T M 1 P  ) 
CS2=ALF0*C0S  CTM2P) 

CS3  =  AL  F 1  *  COS  £ TM  3P ) 
1FCTM.LT.TAU  £1 ) >  C  S 1 =0.0 
IFCTM.LT.TAU £2)  )  CS2=0.0 
IFCTM.LT .TAU £3)  )  CS3  =  0.0 

•  ***  NOISE  GENERATION 

CSN=URAN0CX> 

CSN=FSNR*CSN 
SUM  *S  UM  ♦  C  S  N 
SUM2*SUM2*CSN«*2 
AM  A  G  £  I  )  *  C  S  1 
BMAGC 1 )=CS2 
CMAGC  1 1  * C S  3 
CS*CS1*CS2*CS3 
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YO)  =  CMPLX(CS  .SS) 
AN(I)=CMPLXCC$N,S$N) 
ANS ( 1 ) -CSN 
CONTINUE 


C 

C  *  ***  FINO  MEAN. MEAN  SQUARED, AND  STANDARD  DEVIATION 


FNP*FLOAT (NP ) 
SUM=SUM/FNP 
SUM  2sSUM2/fNP 
STDV=SQRT (SUM 2-SUM* *2) 
C 
c 


c 

C*«**  PRINT  STATISTICAL  PROPERTY  OF  THE  NOISE 
C 


30 


AO 

C 

c 

c 


WRI TE (6  ,2  ) 

WR1TEC6.30) 

FORMAT  (  i  0  X .35HSTAT 1ST ICAL  PROP£RTI£S  OF  THE  NOISE) 
WRITE (6,2) 

WRITE (6, AO)  SUM.SUM2.STDV 

FORMAT (20 X ,5hM£ AN= ,G1 3 .6  /  ,17x ,8HME AN  SQ=,613.6/12X 
♦13HSTANDARD  DEV=,G13.6) 


C*  ***  FIND  F  FT  OF  THE  NOISE  AND  THE  SIGNAL 
C 

CALL  FFTCY.M) 

CALL  FFT(AN.M) 

C 

C 

C“***  FIND  SIGNAL  TO  NOISE  RATIO  AND  PRINT  IT 
C 

TOTLS=0.0 
TOTLN=0.0 
DO  50  1*1  ,NP 

TTS=R  EAL(Y (I ) )  *  *  2  ♦  A  IMAG (Y  (I  )  )  *  *  2 
TTN  =R  E  AL  (  AN  (  I  )  )  *  *  2  ♦  A  I  MA  G  (  AN  (  I  )  )  *  *  2 
TOTLS=TOTLS*TTS 
TOTLN=TOTLN*TTN 
50  CONTINUE 

TOTLS  =  SQRT  (TOTLS) 

TOT  LN  *SQ  R  T (T  0  T  L  N  ) 

SNR=TOTLS/TOTLN 
SNR=10  .*ALOG 1 0 ( SNR  ) 

WRITE (6,60)  SNR 

60  FORMAT (21 X ,4 HSNR* ,F  10 .5 ,5H  DB> 

WRITE (6,2> 

C 


C 

(•«•*  MIX  SIGNAL  WITH  THE  NOISE 
C 

DO  65  1=1, NP 

AM  A  G ( I )=AMAG  (  1  >  *ANS  (1  ) 

DMAG(  I  >=AMAG  (I)*BMAG(I)UMAGCI  ) 
CSX  DM  A  G  (  I  ) 

Y(I  ) *  CMP  L  X (CS,SS) 

65  CONTINUE 
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MP«  300 
WRITE(6,1) 

WR 1 TE (6,70 ) 

70  FORMAT  (10X  ,26HTH1S  IS  GRAPH  OF  THE  NOISE) 

WRITE  (6,2) 

CALL  SHEKL (Q ,AN$ ,  MP  ,H  N A  ME  ,  V NAME  ,  EM  ,  EX , I C) 

WRl TE (6, 1 ) 

WR I TE (6 ,8  0  ) 

80  FORMAT (1QX,33HTH1S  IS  GRAPH  OF  THE  FIRST  SIGNAL) 
WRITE<6,2) 

CALL  SHEKL (Q , AM AG ,MP ,HN AME , VNAM E  ,  EM ,E X , I C) 

WRITE<6,1) 

WRITE (6 ,90) 

90  FORMAT (10X,3AHTH1S  IS  GRAPH  OF  THE  SECOND  SIGNAL) 

WRJ  TE (6  ,2) 

CALL  SHEKL (Q ,BMA6,MP,HNAME, VN AME , E M ,E X , 1 C ) 
IF(TAUC3).EQ.0.0>  GO  TO  105 
WRI  TE  (6 ,1  > 

WRITEC6.100) 

100  FORMAT (10X ,3 3HTH1S  IS  GRAPH  OF  THE  THIRD  SIGNAL) 

WRI TE (6,2) 

CALL  SHEKL (Q ,CMA6»MP,HNAM£,VNAME,EM,£X,1C) 

105  CONTINUE 

WRITE(6»1 ) 

WRITE (6,110) 

110  FORMAT (10X  ,37HTH1S  IS  GRAPH  OF  THE  COMPOSIT  SIGNAL) 
WRITE (6,2) 

CALL  SHEKL (Q ,DMA6,MP,HNAME, VN AME, EM, EX ,IC) 

WRITE  (6,1  ) 

RETURN 

END 
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SUBROUTINE  SIGNAL  (Y) 

COMPLEX  Y (1024) 

COMMON / S IGLOO /F , ALFO , ALF1 , COEF 
COMMON/ S IGCMP / T , NP , TAU ( 3 ) 

C****  THIS  SUBROUTINE  WILL  GENERATE  FEW  DAMPED  COS  SIGNALS 

C****  WITH  THE  SAME  FREQUENCY  BUT  DIFFERENT  STRENGTH 

C****  AND  ARRIVING  AT  DIFFERENT  TIMES 

C****  TAU ( I )  1=1,2,..  IS  THE  DELAY  IN  SECONDS 

C****  F  IS  THE  FREQUENCY  IN  HERTZ 

C 

C 

1  FORMAT  (1H1) 

2  FORMAT  (3 (/) ) 

C 

C****  SET  THE  CONSTANTS 
C 

PI=3, 141592654 

PI2=2  .*PI 

PI2F=PI2*F 

TN=T/NP 

SS=0.0 

CS3=0. 

C 

C 

C 

C****  SIGNAL  GENERATION 
C  DO  20  1=1, NP 

TM=TN*I 
TMI=TM-TAU ( 1 ) 

TM2=TM-TAU (2) 

TM3=TM-TAU (3) 

TMM1=PI2F*TM1 
TMM2=PI2F*TM2 
TMM3=PI2F*TM3 
CSC1=C0S (TMM1) 

CSC2=COS (TMM2) 

CSC3=COS (TMM3) 

CS1=CSC1*EXP(-C0EF*TM1) 

CS2=ALFO*CSC2*EXP(-COEF*TM2) 

CS3=ALF1*CSC3*EXP(-C0EF*TM3) 

IF (TM.LT .TAU (1))CS1=0.0 
IF (TM.LT .TAU (2) ) CS2=0 . 0 
IF (TM.LT. TAU ( 3) )CS 3=0. 

CS=CS1+CS2+CS3 
Y (I) =CMPLX (CS , SS) 

20  CONTINUE 
RETURN 
END 
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SUBROUTINE  SIGNAL(Y) 

COMPLEX  Y  C 10  2  A  ) 

DIMENSION  AMAGC1024),BMAG  (1024)  .CMAGM024) 
DIMENSION  DMAG(1024> 

COMMON/ CM PL2 /NP,TfF,ALF0,ALF1.C0EF,TAU(3) 

COMMON /D AT AA /HNAME ( 4 ) f V NA ME  (4 ) , Q  Cl  0 24  ) 

C 

c***«  THIS  SUBROUTINE  WILL  GENERATE  FEW  DAMPING 
C****  EXPONENTIAL  SIGNALS  STARTING  AT  DIFFERENT  TIMES 
O***  AND  HAVING  DIFFERENT  STRENGTHS 
C 

1  FORMAT { 1  Hi ) 

2  FORMAT(3(/)> 

C 

C  *  ***  SET  THE  CONSTANTS 
C 

SS=0. 

TN=T/NP 
CS3  =0 .0 

c 

c 

c 

C****  SIGNAL  GENERATION 
C 

DO  10  I  =  1  ,  NP 

TM  =  TN  *  I 

TM1 =TM-TAU (1 ) 

TM  2  =TM  -T  AU  C2  > 

TM3=TM-TAU  C3  ) 

TMM1  =  EXP (~COE  F*TMl ) 

TMM2=EXPC-C0EF*TM2) 

TMM3=EXP(-COEF*TM3) 

C  SI =TM1 *TMM1 

CS2-ALF0*TM2*TMM2 

CS3=ALF1 *TM3*TMM3 

1FCTM.LT.TAU (1)  )  CS1=0.0 

1FCTM.LT.TAU  (2))  CS2-0.0 

1FCTM.LT.TAU C3)  )  CS3=0.0 

CS=CS1*CS2*CS3 

AMAGC I )*CS1 

BMAGC I  )=CS2 

CMAGCI  )  =  CS3 

DMAGC I )=CS 

Y(I)=CMPLX(CS»SS) 

10  CONTINUE 
C 

c 

c 

C*  ***  PLOT  THE  SIGNALS 
C 

MP=300 
WRITE  C6 ,1 ) 

WRITE (6.301 

30  FORMAT  (lOX  ,33hTHIS  IS  GRAPH  OF  THE  FIRST  SIGNAL) 
WRITEC6.2) 

CALL  SHEKLCQ , A M AG , MP , HN AM E , VNA M E , 1 . « 1 . *  2 ) 

WRITE (6,1  ) 

WRITE (6,40) 

40  FORMAT  (10X  .34HTHIS  IS  GRAPH  OF  THE  SECOND  SIGNAL) 
WRITE  (6,2) 

CALL  SHEKLCQ  ,BMAG,MP, HNAME, VNAKE,1.*1., 2) 
IF(TAU(3).EQ .0.0)  GO  TO  60 
WRITE(6,1  ) 

WRITE (6.S0) 

50  FORMAT(10X,33HTHIS  IS  GRAPH  OF  THE  THIRD  SIGNAL) 
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WRJTE(6,2)  „  .  . 

CALL  SHEKL(Q ,CMAG,MP,HNAME.VNAME,1  .,1 ..2> 

CONTINUE 
WR1 TE (6,1 > 

FOfiMA7(foxi36H7HIS  IS  GRAPH  Of  THE  COMPDSJT  SIGNAL) 
WRITE (6,2) 

CALL  shekl(q,dmag,mp,hname,vnahe,i.»i.*2> 

WRI TE  (6,1  ) 

RETURN 

END 


\ 

\ 

\ 

\ 
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SUBROUTINE  WIND  OW  (  I  *N  ,  W  ) 

t+  ***  THE  WINOOWING  PROCESS  IS  BASED  ON  KAJ SER-BE SSEL 
C  *  *  * •  (BETA =8.)  EQUATIONS 
S  =  1 .0  E  -8 
1  B  =  N / 2 
A*( I  — I  8) / IB 
X=8.*S0RT  Cl .-A**2 ) 

E  =  1  . 

DE  =  1  •  0 

r*x/2. 

DO  10  J  -1  »  25 
DE=DE* Y/FLOAT ( J) 

SDE=DE**2 
E-E ♦SDE 

IKE*S-SDE)  10,10,20 
10  CONTINUE 

20  X  =  E 

W=X/427.57 

RETURN 

END 
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SUBROUTINE  FFTCX.M) 

COMPLEX  X C1024)  ,u  ,w  ,T 

C  *  •  *  *  THIS  SUBROUTINE  WILL  COMPUTE  DFT  BY  EFT  METHOD 
C****  FOR  REFERENCE  SEE  “THE  FAST  FOURIER  TRANSFORM 
C* *  AND  IT  APPLICATIONS"  BY  JAMES  W.  C OOL E Y  ,  P . A  .  W . L E W I S 
****  AND  PETER  D.  WELCH  IEEE  TRANS  ON  E D »W OL . 1 2  ,  NO . 1 
****  , MARCH  1969  PP  27-34 


1 

2 


3 


4 

5 


N -2 **M 
N2=N/2 
N 1 =N-1 
J=1 

DO  3  1=1  ,N1 
IFCI.GE.J)  60  TO  1 
T*X  (  J  ) 

XCJ)=X(I) 

X ( I ) =T 
K  =N  Z 

JFCK.GE.J)  GO  TO  3 

J-J-K 

K  =  K/2 

60  TO  2 

J=J  *K 

PI=3. 141592653589793 

DO  5  L =1  »  M 

LE=2**L 

LE1=LE/2 

U*.(1.  0,0.0) 

W  =  CMPLX (COS  CP  I /LEI ), SIN  CPI /LEI )> 

DO  5  J=1,LE1 

DO  4  I=J,N,LE 

ID*l*LEl 

T=X(ID)*U 

X(ID)=XCI)-T 

XCI)=X  U)*T 

U  =  U‘W 

RETURN 

END 
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SUBROUTINE  SMOOTH<N,V> 

COMPLEX  YC10Z4) 

THIS  SUBROUTINE  DOES  SMOOTHING  PROCESS  BY  HANNING 
C***»  METHOD 
N 1 =N-1 

X  =  CREALCYCT) )-*REAL  CY(2)))/2  . 
Y(1)=CMPLXCX,A1MAG(YC1>>) 
X*CREAL(Y(NT))>REAL(Y(N>))/2. 
Y(N)=CMPLX(XfAIMAG(YCN))> 

DO  10  1  =  2,  NT 
A  =  REAL  CV  CI-1 )  )/4, 

C  =  REAL  (YC1-»1  )  )  /  4  , 

B*REAL (YCI))/2. 

X=A ♦B ♦ C 

Y(I)=CMPLXCX .AIMAGCYCI))) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  ft  A  NS ET ( M A  X  I N T , NS TR T ) 

COMMON  /MG /RAN  ( 101, GEN(10)tNWRO, BASE, MOD, FBASE.f MOO 
INTEGER  RAN, GEN, BASE, CARRY, REM 

RANSET  AND  URAND  CONSTITUTE  A  MACHINE  INDEPENDENT 
c***»  RANDOM  NUMBER  GENERATOR  hAVING  SPECIALLY  GOOD 
C*««*  STATISTICAL  AND  CORRELATION  P R 0 PE R T I E S , E OUND  IN 
l****  VOLUME  II  Of  THE  FOLLOWING  BASIC  REFERENCE. 

C*  ***  E.G. MCGRATH  .ET  A L , **T E C HN I Q UE S  FOR  EFFICIENT  MONTE 
C»***  CARLO  SIMULATION 

C**t*  VOL  1.  SELECTING  PROBABILITY  DISTRIBUTIONS 
C  *  *  *  *  (AO-762-721) 

C«***  VOL  II.  RANDOM  NUMBER  GENERATION  FOR  SELECTED 
£**•*  PROBABILITY  D  ISTRIBUTIONS (AD-762-722) 

C****  VOL  III.  VARIANCE  REDUCTION  (AD-762-723) 

MAX I=M  AX  INT/4 
18=0 
BAS  E  =  1 

99  IFCfiASE.GT.MAXI)  GO  TO  100 
BAS  E  =  B  AS  £  *4 

IB=  IB*  1 
GO  TO  99 

100  8ASE=2**  IB 
FBA$£=8ASE 
NWft  D  =  4? / I B  *1 

REM  =47-16* (NWRD-1 > 

MO  D  =2  *  *  R  E  M 
FMOO=MOD 
DO  101  N  =  1  *10 
RAN (N)=0 
GEN  CN ) =0 

101  continue 
GENC1 >=5 

DO  200  1=1  ,14 
CAR  R Y  =  0 

DO  190  N  =1 ,NWPD 
GEN(N)=GEN(N)*5*CARRY 
C  AR  R  Y  =  0 

IfCG£N<N).LT .base)  GO  TO  190 
CARRY=GEN (N) /BASE 
GEN(N)=GENCN)-BASE*CARRY 
190  CONTINUE 
200  CONTINUE 

nstart=nstrt 

IFCNSTART.LE .0)  NSTART=2001 
NSTART=2« (NSTART/2)*1 
DO  300  N  =1  ,N  WRD 
NTEMP=NS  tart /base 
R an (N)=N START -N TEMP «BASE 
nstart=ntemp 
300  continue 
RETURN 
END 
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30 


1 


5 


20 


10 


REAL  FUNCTION  U R A  NO  (  F  RAN) 

COMMON /MG /RAN  < 10), G£N(10),NWRD, BASE, MOD, F BASE 
DIMENSION  SUM(IO) 

INTEGER  R AN, GEN,8ASE ,C ARR  Y, SUM  ,PROD,HPROD 

DO  30  1S=1,N*RD 

SUM  C IS )  =  0 

CON  TINUE 

DO  1  I G= 1 ,NW  R D 

N2=NWR D-IG*1 

DO  1  1R=1,N2 

I S= IR ♦ 1G~1 

PROD=RAN(IR)*G£N(IG> 

HPROD=PROD/BASE 
LPROD=PROD-HPROD*BASE 
SUM (IS ) sS  U  M ( IS) ♦LPROD 

1FC1S.LT.NWRD)  SU M (  I S ♦ 1  )  =  SUM C I S ) ♦ HP RO D 

CONTINUE 

N2=NwRD-1 

DO  5  1 S=1 »N2 

CARRY=SUM(IS ) /BASE 

SUM(IS)=SUM(IS) -CARRY*BASE 

SUM  CIS  ♦1  )=SJM(ISM  )  ♦CARRY 

CON  TINUE 

SUM (NWRC)=SUM(NWRD ) -MOD  * ( SUM C NW R D ) / MO D) 

DO  20  IS=1,NWRD 
RAN ( IS )®SUMC IS) 

CONTINUE 
FRAN=SUM (1 ) 

DO  10  IS=2,NWRD 

FRAN=FRAN/FBASE*FLOAT(SUM(IS>> 

CONTINUE 

FRAN  =  FRAN  / FMOD 

URAND= FRAN 

RFTURN 

END 


FMOD 


FUNCTION  CRAN'D 


9904  .8066  .964?  .4108  .5079  .2305  .8254  .6424 


SAMPLE  PROGRAM  FOR  GEf> 
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SUBROUTINE  SH£KL(B,A,NP,XNAME,YNAME,EMl,EXt  *  I C 1 > 
DIMENSION  8(NP),A(NP),LINE(6l),XNAME(4) ,YNAME (A) 

REAL  LINE 

LOGICAL  MANUAL, SYME 

DATA  NI/lHS/  ,L0/1H#  /  ,BLANK/1H  /  ,  D  0  T/ 1  H  .  /  ,  ST  AR  /  1  H*  / 

EMINsEMI 

EMAXsE  XI 

ISC=1C1 

MANUAL=. TRUE . 

SYME=, FALSE. 

IFCEMAX.EQ.EMIN)  manual=. FALSE. 

IFCEMIN.GT.EMAX)  GO  TO  160 
IISC=ISC 

IF(MANUAL)  1ISC=1 
1FCMANUAL)  GO  TO  T 1 
YMIN  =  1  .E38 
YMAXs-1  .£38 
DO  10  1  =  1, NP 

1F(AC1).GT  .YMAX)  YMAX=A(I> 

IFCACD.LT.YM1N)  YMIN=A  Cl  ) 

10  CONTINUE 
GO  TO  19 

11  YMAX=EMAX 
YMIN=EMIN 

19  XMINsABS CBC1 >> 

M  IN  =1 

DO  20  J=2,NP 

IFCA0S  .GE.XM1N)  GO  TO  20 

XM I N- A8S CB  CJ )  ) 

M I N  =  J 

20  CONTINUE 
1FC1ISC.EG.0)  GO  TO  21 
RANGE=YM AX-YMIN 

KAX  1S  =  60.*  C-YMIN/RANGED1  .5 
IF< YMIN.GT  .0  .0)  <AXIS  =  1 
IFC TM AX  .LT .0 .0)  KAX  2  S  a  6 1 
DIS=RANGE/60. 

GO  TO  30 

21  IFC YMIN.GE  .0  .)  GO  TO  22 
IFCYMAX.LE.O.)  GO  TO  23 
KAX  IS  =  31 

SYME=  .TRUE  . 

ABY  =  ABS  C YMIN) 

RANGE=AMAX1CYMAX,ABY) 

D 1 S  =R  ANG  £  /  3G • 

GO  TO  30 

22  KAX  I  S  =  1 
RAN  GE  =  YM  A  X 
OIS=RANGE/60. 

GO  TO  30 

23  KAXIS=6T 
range=-ymin 

DIS=R ANGE/60  . 

30  CONTINUE 

DO  40  K4=1  ,3 
WRITE  (6,50) 

40  CONTINUE 

50  FORMAT(lHO) 

WRI TE (6,60)  XNAME  < 1 ) ,XNAM£  C2) , VNAME  Cl  )  ,YNAME(2)  , 

* YMIN, YMAX 

60  FORMAT (2 X ,A4 ,A4 ,7X ,A4 ,A4 ,3X  ,20HORDINATE  AX1S&  MIN*, 

*G13.6,3x,4hmax*,G13.6) 
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80 


89 

90 


100 

110 

120 


129 

130 


no 

145 

150 


160 

170 


WRITE (6,70)  XNAME (3  )»  XNAME C 4 ) , Y N AM E C 3  )  . V N AM E  C 4 )  ,  D 1 S 

FORMAT  (2  X  ,  A4  ,  A  4  ,  7  X ,  A  4  ,  A  4  , 25X,10HINCREMENT»»613.6/) 

00  80  J  =  1  ,61 

L I N  E ( J ) =8LAN  K 

CONTINUE 

00  100  J2=1,NP 

If< J2.EQ  .MIN)  GO  TO  110 

X=60.0«C(ACJ2)-YMIN)/RAN6E) 

IFCIISC.EQ.O.AND.YMIN.GT.C.)  X=60.*CACj2)/RANGE> 
IFCSYME)  X=30 .* (A ( J2) /RANGE )*30. 

K  =  X  *  1 .5 


LlNE(KAXIS)=DOT 
IFCK.GT.61)  LINE(61)=HI 
IFCK.GT.61)  KMAX=61 
IFCK.GT.61)  GO  TO  89 
IFCK.LT.1 )  l INE ( 1 ) =  LO 
1FCK.LT.1J  KMAX=KAXIS 
IFCK.LT.1 )  GO  TO  89 
LINECK)=STAR 
KM  A  X  =  K  A  X  IS 

1  F C K. GT  .KMAX  )  KMAX=K 

WRIT£(6,90>  B<J2),A(J2),CLINECN4),N4=1,KMAX) 
FORMAT (1X.G13. 6, 2X«G13.6,2X,6lAl) 

IFCK.GT.61)  LINEC61  )=BLANK 

IFCK.GT.61)  GO  TO  100 

IFCK.LT.1)  L INE  (1 ) =8LANK 

IFCK.LT.1  )  GO  TO  100 

LINE(K)=BLANK 

CONTINUE 

GO  TO  145 

CONTINUE 

DO  120  J  = 1  ,61 

LINEC J)=DOT 

CONTINUE 

X=60. * C ( A  CM I N)-YMl N  ) /RANGE) 

IFCIISC.EO.O.ANO.YMIN.GT.O.)  X=60.*CA(J2)/RANGE) 
IFCSYME)  X=30.*(AC J2J/RANGE )*30. 

KP 1 =  X*1  .5 


IFCKP1  .GT.61  )  LINE  (61 ) =H1 
IFCKP1.LT. 1)  LINE  Cl  )  -  LO 
1FCKP1. LT.1. OR. KP1. GT.61)  GO  TO  129 
LINECKPl )=ST AR 

WRITEC6,130)  B(J2),ACJ2),LINE 
FORMAT (lx. G13.6,2X,G13.6,2X,61A1) 

DO  140  J  =  1  ,6  1 
LINEC J)=BLANK 
CONTINUE 
GO  TO  100 
DO  150  K  5  =  1 , 3 
WRI TE (6 ,50) 

CONTINUE 

EMAX=YMAX 

EM1N=YMIN 

RETURN 

CONTINUE 

WRI TE (6 , 1 70) 

FORMAT (10X ,34HHAVE  A  MISTAKE  INPUTTING  EMIN&EMAX) 
RETURN 


F.ND 
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c 


SUBROUTINE  RPRD 

COMMON /S  1GCMP/T,NP#TAU(3),M,SNR 

COMMON  /  CMP LOO/XSZ  ,  YSZ ,  XORG, YORG 

COMMON /CMP LO 1 /Q  CT 0 2 6 )  ,  A  MAG C T 026  )  ,  Y M AX 

COMMON /CMPL02 /I  FOR M.NCOPY  ,LNF  , NOSY M  ,  1 WNDW , I WND F 

GIVE  THE  SIZE  OF  X  AND  Y  NEEDED  FOR  THIS  PLOT 

CALL  PLOTSCXSZ,YSZ  .  IF ORM.NCOPY > 


C 

C  «  «  •  * 

c 

c 

c 

c 

c  *  *  *  * 
c  *  •*  * 

c 

c 

c 

c 

c  •*** 
c 


c 

c 

c 

c*«** 

c 

c 

c 

c 

c 


N=NP 

specify  thickness  of  the  frame 

CALL  LINEWT(-I) 


SET  THE  ORIGIN  AT  THE  POINT  1.,1.  THEN  FROM  NOW 
ON  IT  WOULD  BE  POINT  0.,0. 

CALL  PL0TC1. ,1. ,-3> 


DRAW  THE  FRAME  FOR  THE  SIZE  OF  THE  PAPER 

CALL  PLOT (XS Z  ,0 .0 ,2) 

CALL  PLOT CXSZ,YSZ ,3) 

CALL  PLOT  (0.0  ,  YSZ  ,2  ) 

CALL  PLOT CO.O ,0 .0 ,2 ) 


SET  NEW  ORIGIN 

CALL  PLOT(XORG,YORG,-3) 


C 

c 

c 

c  *  *  *  * 

c 

c 

c 

c 


'  DRAW  THE  COORDINATES  FOR  THE  PLOT 

XCOR=8. 

YCOR  =  5 • 

CALL  SCALECQ Cl) ,XCOR,N,1) 

CALL  SCALE  CAMAG  Cl  >  »YC  OR  ,N  ,1  ) 

N1  =N-»1 
N2=N*2 

AMAG(N2)=YMA  x/ycor 
G(N2)=T/XCOR 

CALL  AXIS(0.0,0.0,I2H  CEPSTRUM  ,  1 2  ,  Y COR  ,9 0 .  , 
•AMAGCM) *A  MA  G ( N2 ) I 

CALL  AX  IS CO.O  ,0  .0 , 1 7HOUEF RE NCY  S E C OND S . -1 7, X C 0 R  , 
*00.0,Q  (N 1  )  t Q CN2  )) 


SPECIFY  THE  THICKNESS  OF  ThE  PLOT 
CALL  LINEwTCO) 


PLOTING  PROCESS 

CALL  HNECQC1  ),AMAG(1  >,NttfLNFtNOSYMfH) 
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C  *  *  *  *  PROCEDURE  TO  WRITE  RELATED  PARAMETERS  ON  THE  PLOT 
C 

CALL  PLOTCO. 0,0. 0,-3) 

W=. 0875 
XW*2.00 

1FCTAU(3).NE  .0.0)  XW  =  1 .5 
VW=-.75 

CALL  SYM80L(XW,YW,H,6H1WNDU*,00.00,6,0) 

XW=XW*6.*W 
FW=FLOAT (IWNDW) 

F  Fs  FLO AT CIWNDF) 

FN=FLOAT  (N  ) 

CALL  NUMBER(XW,YW,H,FW,0.0,-1 ,0) 

XW=XW*W*4  . 

CALL  SYM80L(XW,YW»H,6HIWNDF=, 00.00, 6,0) 

XW*XW*6.*W 

CALL  NUMBER(XW,YW,H,FF,0.0,-1,0> 

XW*XW*W*4. 

CALL  SYM80L(XW,YW,H,2HN=,00  .00  ,2  ,0) 

xw-xu  +  2  •  *  w 

CALL  NUMBER ( XW,YW ,H ,FN ,0. 0,-1 ,0) 

XW=XW*5  .*W 

CALL  SYM80L(XW,YW,H,5HTAUl=, 00.00,5,0) 

XW*  XW  *6  •  *W 

CALL  NUMBER(XW,YW,H,TAU(1  ), 0.0, *3,0) 

XW=XW*6.*W 

CALL  S YMBOLC  XW, YW ,H,5HTAU2=,00 .00 ,5 ,0 ) 

XW*  XW  *6  •  *  W 

CALL  N  UM  BER(XW,YW,H,TAU(2 )#0.Q»*3»0) 

CALL  PLOT CO.,0.  ,-3) 

YCOR=YCOR*1.$*H 

IFCTAUC3) .NE .0.0)  GO  TO  10 

xw=xw*5.*w 

CALL  SYMBOL(XW,YW,H,5H  SN R=  ,00 . 00 , 5  ,0  > 

XW=XW*5.*W  „  „  , 

call  NUMBER£XW,YW,H,SNR,0.0 ,*3,0> 

CALL  SYMBOL(1.0.VCOR,H,74M  C OS ( 2 00 *P I  * ( T -T A Ul  )  > * 

•S(200*PI«(T-TAU2))*UCT-TAU2)*NOISE,0.0,74,0> 

GO  TO  20 
10  CONTINUE 

XW=XW*6.*W 

CALL  SYMBOL(XW,YW,H,5hTAU3=, 00.00, 5,0) 

XW=XW*6.«W 

CALL  NUMBERCXW,VW,H,TAU<3),0.0,*3,0) 

XW*  XW 4 5  . * W 

CALL  S YMBOL( XW, Yw  , H  ,5H  $N  R  =  ,  0  0 . 0  G  ,  5 , 0  ) 

XW*XW*  5  .*  W 


CALL  NUMBER C XU, YW  ,H  ,SNR  .0  .0  ,♦3.0  ) 

CALL  SYMBOLt.5,YCOR,H,l04HCOS(2O0*PI*(T-TAu1>>*U(T 
*0*PI*(T-TAU2))*U(T-TAU2)*.25*COSC200*Pl*(T-TAU3)*U 


*0.0  ,104,0) 

20  CONTINUE 

CALL  PLOT (0.0,0.0,999) 


RETURN 


END 
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KEY  TO  THE  FIGURES 
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This  appendix  presents  time  variation  of  several  signals  which 
have  been  discussed  in  detail  in  the  main  text.  Following  each 
signal  waveform  is  its  cepstral  plot.  In  most  of  the  signal  plots 
or  their  cepstral  plots  the  equation  of  the  composite  signal  is 
written.  If  this  equation  is  broken  in  "+"  points,  each  part 
will  give  the  equation  of  the  individual  waveforms. 

Several  notations  have  been  used  for  the  plots  which  have 
following  meanings. 


IWNOW  =  0  = 
IWNOW  =  1 
IWNOF  =  0  ' 
IWNOF  =  1 
N 

TAl'i  =  ,  . 

l 


no  windowing  process  is  performed  in  time  domain 
windowing  process  is  performed  in  time  domain 
no  smoothing  is  performed  in  frequency  domain 
smoothing  process  if  performed  in  frequency  domain 
total  number  of  sample  points 
beginning  time  for  the  i*"*’  signal 


In  the  following  pages  figures  C- 1  through  G- !i  are  cosine 
waveforms  and  their  cepstral  plots.  The  signals  have  different: 
starting  times  and  also  the  number  of  echoes  which  makeup  the 
somposite  signal  is  different  in  some  eases.  The  equation  lor  the 
signals  is  given  on  each  plot.  In  all  plots  N=1024,  IWN'!)k’=0,  and 
lV.rNI)l'=0  except  if  they  are  defined  otherwise. 
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Figure  0-1  a 
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Figure  G-2a 
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Figure  G-3b .  QUEFRENCY  SECONDS  T  «  .07  T,-  .12 


0,13  0,10  0,25  0.31  0,33  0,44  0.50 

Figure  g-3c .  QUEFRENCY  SECONDS  x  =  .07  T  =  .12 
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Figures  G-6  through  G-10  pictorially  show  damped  cosine  wave¬ 
forms  and  their  cepstral  plots.  The  starting  time  of  the  signals 


as  well  as  the  number  of  them  might  be  different  for  each  case. 

The  equation  for  the  composite  signal  is 

(t-T^)  e  ^  cos  (200tt  (t-Tj)  u  (t-T^) 

+  (t-i2>  e  T2^  cos  (200  n  (t-i,,))  u (t-i2) 

+  (t-T^)e  T3^  cos  (200  ti  (t-T^))  u  (t-T^) 

with  the  third  component  being  zero  where  only  one  echo  signal  has 
been  studied.  If  the  above  equation  is  broken  in  "+"  points  each 
part  will  give  the  equation  of  the  individual  signals.  In  all 
figures,  1WNDW=~,  IWNDF=0,  and  N=1024  except  if  they  are  defined 
otherwise . 
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Figure  G-8e .  OUEFRENC'I  SECONDS  T  =  .07  T  -  .12 


TIME  SECONDS 
0BMPEO  COSINE 


T2  =  .09  Sec 


orupeo  co5x„e 

Figure  C-]0,i 


Figures  G-ll  through  G-15  present  the  damped  exponential 
waveforms  with  their  cepstral  plots.  The  damping  coefficient  for 
all  of  these  plots  is  30  ,  and  the  whole  period  of  study  is  .5 
second.  The  equation  of  the  composite  signal  is  given  in  each 
plot,  which  if  broken  in  "+"  points  will  result  in  giving  the 
equation  of  each  individual  signal.  In  all  figures  N=1024, 
IWNDW=0,  and  IWNDF=0  except  if  they  are  defined  otherwise. 
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Figure  r,-15b  qUEFRENCY  SECONDS  t  =  .07  x  =  .09  t  =  .15 
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In  the  following  pages  the  following  figures  have  different 
significances.  Figure  G-16  shows  the  effect  of  under  sampling,  its 
sampling  rate  is  452  samples/sec.  It  does  not  represent  the  cosine 
waveform  in  its  exact  shape.  Figure  G-16  shows  the  cepstrum  plot 
of  Figure  G-16a,  which  seems  to  be  noisy.  Figure  G-16e  is  the 
enhanced  version  of  Figure  G-16b,  meaning  that  after  taking  the 
Fourier  transform  of  the  cosine  signal  for  the  first  time,  the 
value  of  two  impulses  have  been  set  equal  to  one.  As  we  see 
from  Figure  G-16c  there  is  some  improvement. 

Figures  G-17a  and  G-17b  show  effect  of  having  different 
strength  for  the  echoes.  As  these  Figures  show  there  is  no 
difference  in  the  cepstral  property  of  the  plots  except  that  the 
strength  of  the  peaks  has  been  changed  accordingly. 

Figures  G-18a,  G-18b,  and  G-18c  are  the  cepstral  plots  when 
there  are  three  echoes  besides  the  original  signal.  As  is  obvious 
from  these  cepstral  plots,  it  will  be  very  hard  to  distinguish  the 
exact  locations  of  each  epoch  or  echo  delay  time.  In  all  figures 
N=1024,  IWNDW=0,  and  IWNDF=0  except  if  they  are  defined  otherwise. 


.cnsi2oo»p 


Figure 


181 


m 

H 


- *_ 

e*'on 


wnyisdao 


AD-A102  358  MISSISSIPPI  STATE  UNIV  MISSISSIPPI  STATE  ENGINEER I M6— ETC  F/G  17/1 

INVESTIGATION  OF  CEPSTRUM  ANALYSIS  FOR  SEISMIC/ACOUSTIC  SIGNAL  — ETC(U) 
JAN  81  F  M  INGELS »  G  KOLEYNI  AFOSR-80-OOB6 

UNCLASSIFIED  MSSU-EIRS-EE-8l-a  AFOSR-TR-81-0603  NL 


EXP0NENTIRL 


189 


unuisd33 


L 


sr  '  1  60'  ■  *  l0'  ■  1  S0N033S  X3N3U  J3(l6  ',0M  ai'’3,J 

0S*0  tvo  SE*Q  1£*0  £E*0  61*0  £1*0  50*0  CO*(L 

l _  __ _ _ ( _ _ _ _ _ f _ _ _ _ I _ _ i _ i _ _  _ _ _ . _ _ ■ _ i  ® 


APPENDIX  H 


BIBLIOGRAPHY  AND  LITERATURE  SURVEY 


LITERATURE  SURVEY 


In  an  effort  to  appraise  the  state  of  development  of  cepstrum 
technique,  an  extensive  survey  of  cepstrum  related  literature  was 
undertaken.  The  dates  of  the  survey  extends  from  the  beginning  of 
publication  of  a  paper  by  Bogert,  et  al,  in  1963  to  the  present 
time.  The  majority  of  this  literature  was  available  in  the  library 
of  the  Mississippi  State  University,  or  have  been  furnished  to  the 
author  through  the  Inter  Library  Loan.  The  sources  of  the  articles 
and  the  inclusive  search  dates  are  as  follows: 


1.  Dissertation  Abstract  International,  B,  The  Science  and 
Engineering,  January  1963  through  August  1980. 

2.  Geophysics,  January  1960  through  August  1980. 

3.  Geophysics  Journal  of  Royal  Astronomical  Society,  1970 

A.  Institute  of  Electrical  and  Electronics  Engineers  Proceedings, 
January  1963  through  August  1980. 

5.  Institute  of  Electrical  and  Electronics  Engineers  Spectrum, 
January  1963  through  August  1980. 

6.  Institute  of  Electrical  and  Electronics  Engineers  Transactions 
on  Acoustic  Speech  and  Signal  Processing  formerly  Institute 

of  Electrical  and  Electronics  Engineers  Transactions  on  Audio 
and  Electroacoustics,  January  1963  through  August  1980. 

7 .  Institute  of  Electrical  and  Electronics  Engineers  Transactions 
On  Communications,  January  1963  through  August  1980. 

8 .  Institute  of  Electrical  and  Electronics  Engineers  Transactions 
on  Education,  March  1969. 

9.  Institute  of  Electrical  and  Electronics  Engineers  Transactions 
on  Information  Theory,  January  1963  through  August  1980. 

1 0 .  Institute  of  Electrical  and  Electronics  Engineers  Transactions 
on  Systems  Man  and  Cybernetics,  January  1972. 


i. 
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11.  Journal  of  Acoustical  Society  of  America,  January  1963 
through  August  1980. 

12.  Journal  of  Geophysical  Research,  October  1970. 

13.  Journal  of  Sound  and  Vibration,  January  1974  through  August 
1980. 

14.  Report  on  8th  Annual  Southeastern  Symposium  on  System  Theory, 
April  1976. 

15.  Simulation,  March  1969. 

16.  The  Bell  System  Technical  Journal,  January  1963  through 
August  1980. 

For  the  purpose  of  making  a  meaningful  summary,  the  afore¬ 
mentioned  topics  will  be  broken  in  two  categories:  (1)  Fourier 
transform  ,  specially  the  Fast  Fourier  Transform  (FFT) ,  (2)  Cepstrum 

or  the  related  topics.  Included  in  the  summary  of  each  category 
will  be  a  list  of  the  most  pertinent  articles  with  a  synopsis. 
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FOURIER  TRANSFORM 

The  fourier  transform  of  f (x)  is  defined  as 

CO 

r  r /  S  -i 2wxs  , 

J  f  (x)  e  J  dx 

—  CO 

this  integral,  which  is  a  function  of  s,  may  be  written  F(s). 
Transforming  F(s)  by  the  same  formula,  we  have 

CO 

r  .  -j2irws, 

J  F(s)  e  ds 

— oo 

when  f (x)  is  an  even  function  of  x,  that  is,  when  f(x)  =  f (-x) , 
the  repeated  transformation  yields  f (w) ,  the  same  function  with 
which  we  began.  This  is  a  cyclical  property  [82]  of  the  Fourier 
transformation,  and  since  the  cycle  is  of  two  steps,  the  reciprocal 
property  is  implied:  if  F(s)  is  the  Fourier  transform  of  f (x) , 
then  f(x)  is  the  Fourier  transform  of  F(s). 

The  cyclical  and  reciprocal  properties  are  imperfect,  however, 
because  when  f(x)  is  odd  i.e.,  when  f (x)  =  -f(-x)'  the  repeated 
transformation  yields  f(-w).  In  general,  whether  f(x)  is  even  or 
odd  or  neither,  repeated  transformation  yields  f(-w). 

The  customary  formula  exhibiting  the  reversibility  of  the 
Fourier  transformation  are  [82] 

CXI 

F(s)  =  /  f(x)  e-j2nsdx 

—  CO 

f(x)  =  J  F(S)  ej2,‘SX  ds  . 

In  this  form,  two  succesive  transformations  are  made  to  yield  the 


original  function. 
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It  may  happen  that  function  values  are  given  only  at  discrete 
values  of  the  independent  variable  (pulses) ,  as  which  physical 
measurements  made  at  regular  time  intervals.  Regardless  of  the 
form  of  the  given  function,  if  the  transform  is  evaluated  by 
numerical  computing,  the  values  of  the  transform  will  be  available 
only  at  discrete  intervals.  We  often  think  of  this  as  though  an 
underlying  function  of  a  continuous  variable  really  exists  and  we 
are  approximating  it. 

The  rules  for  finding  the  frequency  response  samples  from 
the  pulse  response  and  vice-versa  were  called  the  discrete  Fourier 
Transform  (DFT)  and  Inverse  Discrete  Fourier  Transform  (IDFT) , 
respectively  [81]. 

The  Discrete  Fourier  Transform  or  DFT  of  the  N-point  sequence 

f  ,f  ,...f  is  defined  as  the  N-point  sequence. 

U  1  1 


N-l 

I 


n=0 


f 

n 


-j2"n 


k 

N 


The  inverse  discrete  Fourier  transform  formula  or  IDFT  is 

f  =  IDFT  [F  ]  =  ^  \  Fkej2'Tn-^- 

n  k  N  ,  n  k  N 

k=u 


n  =  0,1,2, ...N-l 


In  1965  a  method  of  computing  discrete  Fourier  transforms 
suddenly  became  widely  known  (J.W.  Cooley  and  J.  W.  Tukey,  math, 
comput,  vo 1 .  19,  April  1965,  pp.  297-301)  which  revolutionized  many 
fields  where  onerous  compting  was  an  impediment  to  progress. 
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Following  is  a  chronological  list  of  some  of  the  most 
pertinent  articles  related  to  Fourier  transform  (FT)  and  Fast 
Fourier  Transform  (FFT) .  A  brief  synopsis  is  included  with  each 
article  in  this  list. 

1.  Cochran,  W.  T.,  Cooley,  J.  W. ,  Favin,  D.  L. ,  Helms,  H.  D. , 
Kaenel,  R.  A.,  Lang,  W.  W. ,  Maling,  Jr.  0.  C.,  Nelson,  D.  E. , 
Rader,  C.  M.  and  P.  W.  Welch,  "What  is  Fast  Fourier  Transform?", 
IEEE  Transactions  on  Audio  and  Electroacoustics,  Vol.  Au-15, 

N8.  2,  pp.  45-55,  June  1967. 

In  this  paper,  the  discrete  Fourier  transform  of  a 
time  series  is  defined,  some  of  its  properties  are  dis¬ 
cussed,  the  associated  fast  method  (Fast  Fourier  Transform) 
for  computing  this  transform  is  derived,  and  some  of  the 
computational  aspects  of  the  methods  are  presented. 

2.  Cooley,  J.  W.  ,  Lewis,  P.  A.  W.  and  P.  D.  Welch,  "Historical 
Notes  on  the  Fast  Fourier  Transform",  IEEE  Transactions  on 
Audio  and  Electroacoustics,  Vol.  Au-15,  No. 2,  pp. 76-84, 

June  1967. 

In  this  paper,  the  contributions  of  many  investigators 
are  described  and  placed  in  historical  perspective. 

3.  Rader,  C.  M. ,  "Discrete  Fourier  Transforms  When  the  Number  of 
Data  Samples  is  Prime",  IEEE  Proceedings,  Vol.  56,  No.  5, 

pp.  1107-1108,  June  1968. 

The  discrete  Fourier  transform  of  a  sequence  of  N  points 
where  N  is  a  prime  number,  is  shown  to  be  essentially  a 
circular  correlation.  This  can  be  recognized  by  rearranging 
the  members  of  the  sequence  and  the  transform  accoring  to  a 
rule  involving  a  primitive  root  of  N.  This  observation  permits 
this  discrete  Fourier  transform  to  be  computed  by  means  of  a 
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fast  Fourier  transform  algorithm,  with  the  associated 
increase  in  speed,  even  through  N  is  prime. 

4.  Cooley,  J.  W. ,  Lewis,  P.A.V.  and  P.  D.  Welch,  "The  Fast  Fourier 
Transform  and  its  Applications",  IEEE  Transactions  on  Education, 
Vol .  12,  No.  1,  pp.  27-34,  March  1969. 

A  description  of  the  algorithm  for  FFT  and  its 
programming  is  given  followed  by  a  theorem  relating 
its  operands,  the  finite  sample  sequences,  to  the 
continuous  functions  they  often  are  intuned  to 
approximate.  An  analysis  of  the  error  due  to 
discrete  sampling  over  finite  ranges  is  given  in 
terms  of  aliasing.  Procedures  for  computing  Fourier 
integrals,  convolution  and  lagged  products  are  outlined. 

5.  Bergland,  G.  D.,  "A  Radix-Eight  Fast  Fourier  Transform 
Subroutine  for  Real-Valued  Series",  IEEE  Transactions 
on  Audio  and  Electroacoustics ,  Vol.  AU-17,  No.  2, 

pp.  138-144,  June  1969. 

Fast  Fourier  Analysis  (FFA)  and  fast  Fourier 
synthesis  (FFS)  algorithms  are  developed  for  computing 
the  discrete  Fourier  transform  of  a  real  series,  and 
for  synthesizing  a  real  series  from  its  complex  Fourier 
coefficients. 

6.  Singleton,  R.  C.,  "A  Short  Bibliography  on  the  Fast  Fourier 
Transform” ,  IEEE  Transactions  on  Audio  and  Electroacoustics , 

Vol.  Au-17,  Mo.  2,  pp.  1 66- 1 69, Jan.  1969. 

A  chronologically  listed  name  of  papers  on  Fast  Fourier 
Transform",  nn.  166-169,  January  1C,69. 
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7.  Singleton,  R.  C. ,  "An  Algorithm  for  Computing  the  Mixed  Radix 
Fast  Fourier  Transform",  IEEE  Transactions  on  Audio  and  Elec¬ 
troacoustics  ,  Vol.  Au-17,  No.  2,  pp.  93-103,  June  1969. 

The  paper  presents  an  algorithm  for  computing  the 
fast  Fourier  transform  based  on  a  method  proposed  by  Cooley 
and  Tukey.  As  in  their  algorithm,  the  dimension  n  of 
transform  is  factored  (if  possible),  and  n/p  elementary 
transforms  of  dimension  are  computed  for  each  factor  p 
of  n.  The  algorithm  is  described  mathematically  and 
illustrated  by  a  Fortran  subroutine. 

8.  Uhrich,  M.  L. ,  "Fast  Fourier  Transforms  Without  Sorting", 

IEEE  Transactions  on  Audio  and  Electroacoustics,  Vol.  Au-17, 

No.  2,  pp.  170-172,  June  1969. 

9.  Markel,  J.  D. ,  "FFT  Pruning”,  IEEE  Transactions  on  Audio  and 
Electroacoustics ,  Vol.  Au-19,  No.  4,  pp.  305-311,  December  1971. 

There  are  basically  four  modifications  of  the  N=2m  point 
FFT  algorithm  developed  by  Cooley  and  Tukey  which  give 
improved  computational  efficiency.  It  is  shown  that  for 
situations  in  which  the  relative  number  of  zero-valued 
samples  is  quite  large,  significant  time  saving  can  be 
obtained  by  pruning  the  FFT  algorithm. 

10.  Harris,  F.  J.,  "On  the  Use  of  Windows  for  Harmonic  Analysis 
With  the  Discrete  Fourier  Transform",  TEEE  Proceedings,  Vol.  66, 
No.  1,  pp.  51-83,  January  1978. 

This  paper  makes  available  a  concise  review  of  data 
windows  and  their  effect  on  detection  of  harmonic  signals 


in  the  presence  of  broad-band  noise,  and  in  the  presence 
of  nearby  strong  harmonic  interference.  Also  calls  attention 
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to  a  number  of  common  errors  in  the  application  of  windows 
when  used  with  the  Fast  Fourier  Transform. 

In  addition  to  the  ten  articles  listed  above,  the  following 
references  to  the  bibliography  in  Appendix  A  are  also  related  to 
Fourier  transform  or  Fast  Fourier  Transform  (FFT) :  1,12,15,16,61, 
66,69, 70, 77, 78, 79, 80, 81, and  82. 


CEPSTRUM  AND  DECONVOLUTION 


The  first  paper  on  the  cepstrum  has  been  published  by  Bogert, 
et  al.  [3]  in  1963  and  this  has  been  the  first  paper  in  inventing 
and  using  this  word  a  thorough  investigation  of  the  subject  has 
been  done  through  the  text.  The  papers  mainly  deal  with  problems 
of  power  cepstrum.  Also  included  are  papers  related  to  the  complex 
cepstrum  as  defined  by  Schaffer  [16].  The  following  list  of 
articles  applies  to  cepstrum  analysis  and  related  subjects. 

1.  Bogert,  B.  P.,  Healy,  M.  J.  R.  and  J.  W.  Tukey,  "The 
Quefrency  Analysis  of  Time  Series  for  Echoes:  Cepstrum, 
Pseudo-Autocovariance,  Cross-Cepstrum  and  Saphe  Cracking", 
in  Pro.  Symp.  on  Time  Series  Analysis,  M.  Rosenblatt,  Ed., 

New  York,  Wiley,  Chap.  15,  pp.  209-243,  1963. 

The  first  published  paper  on  cepstrum  (power  cepstrum) 
brings  up  the  question  of  echoes  and  how  to  detect  them,  also 
compares  the  cepstrum  method  with  autocorrelation.  The 
authors  define  new  terms  paraphrased  from  known  terminologies 
such  as  cepstrum  which  is  a  paraphrased  word  from  spectrum. 
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2.  Noll,  A.  M. ,  "Short-Time  Spectrum  and  Cepstrum,  Techniques 
for  Vocal-Pitch  Detection",  J.  Acoust.  Soc.  Am.,  Vol.  36, 

No.  2,  pp.  296-302,  February  1964. 

In  this  paper  the  author  is  trying  to  use  the  newly 
familiar  subject  of  cepstrum  for  vocal-pitch  detection  in 
speech.  The  cepstrum  of  a  speech  signal  has  a  peak 
corresponding  to  the  fundamental  period  for  voice  speech 
but  no  peak  for  unvoiced  speech.  Thus  a  cepstrum  analyzer 
can  function  both  as  a  pitch  and  as  voiced-unvoiced 
detector. 

3.  Noll,  A.  M. ,  "Clipstrum  Pitch  Determination",  J.  Acoust. 
Soc.  Am.,  Vol.  44,  No.  6,  pp.  1585-1591,  December  1968. 

In  this  paper  the  author  presents  a  new  method  of 
pitch  Determination  similar  to  the  cepstrum  except  that 
both  the  time  signal  and  the  log  power  spectrum  are 
infinitely  peak  clipped  before  spectrum  analysis  has 
been  simulated  on  a  digital  computer.  It  shows  the 
clipping  in  the  clipstrum  might  offer  some  advantages 
over  cepstrum  analysis  in  certain  digital  hardware 
implementations  since  multiplications  could  be  replaced 
with  addition  or  subtractions. 

4.  Childers,  D.  C.,  Varga,  R.  S.  and  N.  W.  Perry,  Jr., 
"Composite  Signal  Decomposition",  IF.EF,  Transactions  on 
Audio  and  F.lect roacoust ics ,  Vol.  Au-18,  No.  4,  pp.  461-477, 
December  1970. 

This  paper  presents  a  technique  for  decomposing 
a  composite  signal,  which  consists  of  the  superposition 
of  known  multiple  signals  overlapping  in  time,  digital 
Data  processing  problems  such  as  filter  realizability. 


signal  resolution  capability,  the  effect  of  additive 
noise,  frequency  (spectrum)  compatibility  between 
signal  waveform  and  filter  response  pulse,  and  possible 
additional  processing  in  certain  cases  are  discussed. 

Kemerait,  R.  C. ,  and  D.  G.  Childers,  "Signal  Detection  and 
Extraction  by  Cepstrum  Techniques",  IEEE  Transactions  on 
Information  Theory,  Vol.  IT-18,  No.  6,  pp.  745-759,  November 
1972. 

Digital  data-processing  problems  such  as  the 
detection  of  multiple  echoes,  various  methods  of  linear 
filtering  the  complex  cepstrum  the  picket-fence  phenomenon, 
minimum-maximum  phase  situation,  and  amplitude-versus 
phase-smoothing  for  additive-noise  cases  are  examined 
empirically  and  where  possible  theoretically,  and  are 
discussed. 

Hassab,  J.  C. ,  "Time  Delay  Processing  Near  the  Ocean  Surface' 
J.  Sound  and  Vibration,  Vol.  35,  No.  4,  pp.  489-501,  1974. 

A  study  of  a  simple  multi-path  channel  near  the  ocean 

surface  where  a  direct  path  and  a  surface-reflected  path 

link  the  source  to  the  receiver.  By  using  autocorrelation 

and  cepstrum  processing,  extraction  of  the  difference  in 

arrival  times  between  two  paths,  i.e.  time  delay,  from  a 

priori  unknown  signals  is  studied. 

Hassab,  J.  C.  and  R.  Boucher,  "A  Probabilistic  Analysis  of 
Time  Delay  Extraction  by  Cepstrum  in  Stationary  Gaussian 
Noise",  IEEE  Transactions  on  Information  Theory,  Vol.  IT-22, 
No.  4,  July  1976. 
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A  probabilistic  analysis  is  conducted  dealing  with 
the  effect  of  stationary  Gaussian  noise  on  the  characteristics 
of  such  a  nonlinear  processor.  The  expected  mean  and 
standard  deviation  of  reduction  in  the  peak  level  at 
(time  delay)  due  to  noise  are  analytically  described. 

The  results  point  out  the  dependence  of  statistical 
measures  upon  the  pointwise  variation  of  input  signal 
to  noise  spectra. 

8.  Childers,  D.  G. ,  Skinner,  D.  P.  and  R.  C.  Kemerait,  "The 
Cepstrum:  A  Guide  to  Processing",  IEEE  Proceedings,  Vol.  65, 
No.  10,  October  1977. 

This  paper  is  a  pragmatic  tutorial  review  of  the 
cepstrum  literature  focusing  on  data  processing.  The 
effects  of  various  forms  of  liftering  the  cepstrum  are 
described.  The  results  obtained  by  applying  witening 
and  trend  removal  techniques  to  the  spectrum  prior  to 
the  calculation  of  the  cepstrum  are  discussed.  A  good 
many  numbers  of  references  are  given  in  this  paper. 

9.  Dudgeon,  D.  E. ,  "The  Computation  of  Two-Dimensional  Cepstra", 
IEEE  Transactions  on  Acoustics,  Speech, and  Signal  Processing , 
Vol.  ASSP-25 ,  No.  6,  December  1977. 

This  paper  explores  two  methods  of  computing  the 
complex  cepstrum  of  a  two-dimensional  signal.  It  considers 
the  definitions  of  two-dimensional  causality  and  two- 
dimensional  minimum  phase  signals.  It  explores  the  relation¬ 
ship  among  the  nonzero  regions  of  a  signal,  its  inverse,  and 


Its  cepstrum. 
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10.  Ingels,  F.  M. ,  "Cepstrum  Analysis  Techniques  for  Possible 
Application  to  Seismic/Acoustic  Ranging",  Report  No. 
F‘-'>9620-79-C-0038 ,  Mississippi  State  University,  August  1979. 

This  report  is  concerned  with  ranging  problems.  It 
studies  use  of  one  sensor  to  determine  the  range.  It 
advocates  cepstral  analysis  to  determine  the  differences 
in  time  of  arrival  of  two  signals  detected  with  the  same 
sensor,  also  shows  some  results  of  applying  the  cepstrum 
technique. 

Other  papers,  articles,  letters,  or  books  which  are  related 
to  the  cepstrum  technique  are  listed  in  the  bibliography  of 
Appendix  A.  The  following  list  of  numbers  indicate  those  entries 
in  the  bibliography:  1,5,6,8,9,10,11,12,16,25,26,27,29,31,32,34, 
35,37,39,40,42,43,44,46,47,48,49,50,51,52,53,54,55,56,57,58,60, 
67,71,73,74,77,79,80,  and  81. 

A  COMPILED  BIBLIOGRAPHY  ON  CEPSTRUM  AND  RELATED  TOPICS 

The  following  bibliography  concerns  the  subjects  of  cepstrum 
homomorphic  deconvolution,  Fourier  transform.  Fast  Fourier  Transform 
(FFT),  power  spectrum  and  applications.  The  bibliography  is 
divided  into  two  sections  as  follows:  (1)  articles  published  in 
journals,  (2)  books.  The  entries  are  listed  chronologically 
except  when  the  month  and  the  year  of  some  publications  are  the 
same  in  which  case  the  listing  is  alphabetical.  An  alphabetical 
author's  index  follows  the  publications  listing  with  reference  to 


numbered  articles  of  the  bibliography. 
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